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Abstract 

We  have  collected  and  analyzed  large  and  diverse  datasets  of  seismic  waveforms,  absolute 
and  relative  body  wave  travel  times,  and  surface  wave  phase  velocities  and  amplitudes. 
Our  primary  objective  has  been  to  map  the  variations  of  elastic  mantle  properties  beneath 
Eurasia  over  horizontal  lengthscales  of  approximately  1000-1500  km  and  vertical  lengthscales 
of  50-100  km  in  the  upper  mantle. 

In  order  to  characterize  and  properly  account  for  the  highly  heterogeneous  velocity  struc¬ 
ture  associated  with  the  crust  and  uppermost  mantle,  we  have  made  several  tens  of  thousands 
of  measurements  of  intermediate  period  surface  wave  dispersion,  and  have  developed  global 
maps  of  phase  velocity  for  Love  and  Rayleigh  waves  in  the  period  range  35-150  s.  The  new 
maps,  expanded  in  spherical  harmonics  up  to  degree  40,  provide  up  to  96%  reduction  of  the 
variance  observed  in  the  phase  measurements,  and.  at  short  periods,  show  a  very  strong  cor¬ 
relation  with  the  thickness  of  the  crust  and  the  surface  geology.  The  new  dataset  of  surface 
wave  dispersion,  which  provides  excellent  constraints  on  crustal  and  lithospheric  structure, 
has  been  combined  with  previously  collected  data  in  the  finalization  of  our  new  3-D  whole 
mantle  model.  In  comparison  with  previous  Harvard  mantle  models,  such  as  S12WM13,  we 
have  increased  the  horizontal  resolution  by  expanding  the  global  structure  up  to  spherical 
harmonic  degree  I  =  20,  and  we  have  adopted  a  radial  parameterization  which  gives  better 
vertical  resolution  in  the  upper  mantle. 

The  new  model  shows  many  of  the  same  features  seen  in  previous  tomographic  images,  but 
with  significantly  sharper  definition,  in  particular  in  the  upper  mantle.  Lithospheric  roots 
beneath  the  stable  continental  interiors  are  the  most  prominent  feature  in  the  uppermost 
mantle.  A  surprising  and  complicating  result  of  our  study  is  that  when  P/SV-  and  SH- 
sensitive  data  are  inverted  separately,  differences  in  the  resulting  models  are  apparent  in 
the  middle  of  the  upper  mantle  (200-300  km).  The  suggestion  is  that  anisotropy  may  be  a 
significant  contributor  to  lateral  heterogeneity,  particularly  beneath  the  oceans. 

We  have  made  efforts  to  validate  our  models,  and  to  explore  their  utility  in  event  loca¬ 
tion,  through  the  prediction  of  Pn  travel  times.  In  some  regions  we  obtain  relatively  good 
agreement  between  the  observed  residual  and  the  Pn  travel  time  predicted  from  the  P  and 
S  model  S&P12WM13  (expanded  only  to  1=12).  In  other  regions  the  agreement  is  poor, 
probably  indicating  limited  resolution  of  the  structure  in  the  uppermost  mantle. 
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Objective 

We  adapt  the  seismic  tomographic  techniques  developed  in  global  scale  studies  to  the 
regional  scale  problem  of  mapping  the  elastic  and  anelastic  material  properties  beneath 
Eurasia.  The  main  objective  of  this  research  is  to  obtain  a  tomographic  seismic  velocity 
model  for  the  structure  beneath  Eurasia  with  a  horizontal  resolution  corresponding  to  at 
least  I  =  20,  consistent  with  the  widest  possible  range  of  seismological  observations. 

Use  of  tomographic  models  of  the  mantle  in  the  calculation  of  teleseismic  travel  times  have 
been  shown  to  result  in  improved  teleseismic  event  locations  {Smith  and  Ekstrom,  1996). 
These  models  are  therefore  relevant  to  the  global  monitoring  of  a  Comprehensive  Test  Ban 
Treaty,  in  particular  for  obtaining  better  locations  for  events  in  previously  uncalibrated 
regions.  In  addition,  the  improved  ability  to  model  dispersed  intermediate-period  surface 
waves  makes  it  possible  to  perform  more  efficient  filtering  and  stacking  and,  as  a  consequence, 
to  achieve  both  a  lower  surface  wave  detection  threshold,  and  to  obtain  better  surface  wave 
magnitudes.  The  advance  in  our  ability  to  model  surface  waves  can  therefore  benefit  both 
the  detection  and  discrimination  tasks  of  a  monitoring  program. 

Research  Accomplished 

In  earlier  reports  (1993,  1994,  1995)  we  have  presented  preliminary  results  of  our  higher 
resolution  upper  mantle  model  S20U7L5.  We  have  also  shown  the  improvement  of  teleseismic 
event  locations  using  heterogeneous  mantle  models.  In  this  report  we  discuss  (1)  some  results 
for  the  surface  wave  dispersion  characterization  at  intermediate  periods;  (2)  new  indications 
for  anisotropy  in  the  upper  mantle  which  have  become  apparent  as  we  finalize  the  model 
S20U7L5.  and;  (3)  modeling  of  Pn  residuals  using  global  tomographic  models.  An  extensive 
report  on  the  surface  wave  dispersion  results  (1).  in  the  form  of  a  preprint,  is  included  as  an 
appendix. 

Surface  Wave  Dispersion 

Due  to  their  sensitivity  to  the  elastic  structure  of  the  crust  and  the  top  of  the  mantle, 
measurements  of  intermediate-period  surface  wave  dispersion  are  ideal  for  imaging  the  shal¬ 
lowest  portions  of  the  Earth.  The  main  difficulty  in  obtaining  dispersion  measurements  from 
teleseismically  observed  surface  waves  at  periods  less  than  about  100  s  is  that  variations  in 
phase  velocities  are  sufficiently  big  to  cause  path  variations  of  several  cycles  in  phase,  and 
isolated  residual  phase  measurements  [— tt.tt]  at  a  single  period  are  therefore  ambiguous  in 
multiples  of  2~.  For  minor  arc  observations,  phase  measurements  at  periods  longer  than 
100  s  can.  however,  be  associated  with  a  total  cycle  count.  If  a  continuous  dispersion  curve 
can  be  anchored  at  long  periods,  the  total  plmse  perturbation  can  then  be  inferred  without 
ambiguity. 

In  order  to  obtain  the  complete  broadband  dispersion,  we  have  developed  a  method  we 
call  iterative  frequency  band  expansion.  We  initially  determine  the  dispersion  parameters  in 
a  narrow  frequency  range  centered  around  100  s,  and  then,  in  several  iterations,  increase 
the  higher  frequency  end  of  the  total  frequency  range.  In  each  iteration,  the  dispersion 
parameters  are  adjusted  to  minimize  the  residual  dispersion  between  an  observed  and  a 
model  surface  wave  using  the  downhill  simplex  method  of  N elder  and  Mead  (1965).  Once 
a  misfit  between  the  two  wavegroups  is  minimized  for  a  given  frequency  range,  the  range  is 
expanded  to  include  higher  frequencies.  Because  of  the  smoothness  of  the  dispersion  curve,  a 
prediction  can  be  made  for  frequencies  slightly  higher  than  the  ones  included  in  the  previous 
iteration.  This  ensures  that  the  phase  difference  between  model  seismogram  and  observed 
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seismogram  always  remains  small,  so  that  gradient  methods  continue  to  work  in  the  next 
iteration. 

This  method  of  analysis  has  been  systematically  applied  to  digital  seismograms  from  the 
IRIS/USGS,  IRIS/IDA,  GEOSCOPE,  CDSN,  GTSN,  and  MEDNET  networks,  using  earth¬ 
quake  sources  for  the  period  January  1989-September  1995.  Analysis  was  attempted  for 
all  earthquakes  during  this  time  period  with  M\v  >  5.5.  We  use  the  moment  tensors  and 
centroid  locations  in  the  Harvard  GMT  catalog  to  calculate  the  source  excitation.  Only  shal¬ 
low  earthquakes  {h  <  50  km)  were  included,  in  order  to  maximize  the  excitation  of  surface 
waves.  In  addition,  we  only  consider  paths  in  the  distance  range  25°  <  A  <  150°,  in  order 
to  avoid  problems  of  separation  of  the  fundamental  mode  at  short  distances  and  close  to  the 
antipode. 


Figure  1:  Map  shows  the  phase  velocity  variations  of  35  s  Love  waves  for  Eurasia.  Note  the 
strong  correlation  between  slow  velocitit's  and  thick  and  low- velocity  crust. 

More  than  128.000  station-receiver  paths  were  analyzed.  Approximately  37,000  high- 
quality  meiisurements  for  Rayleigh  waves  at  50  s  resulted  from  the  analysis,  while  at  35  s, 
the  number  was  somewhat  smaller  —  28.000.  Fewer  good  measurements  were  obtained  for 
Love  waves,  in  part  due  to  the  higher  noise  levels  on  horizontal  components.  In  the  analysis 
of  the  mea.surements,  the  dispersion  curves  were  converted  to  phase  anomalies  5$  at  the 
discrete  periods  150,  100,  75,  60,  50.  45.  40.  37,  and  35  s.  Observational  uncertainties  of 
the  measurements  were  determined  from  statistical  analysis  of  repeated  measurements  for 
similar  paths. 

The  phase  anomalies  were  next  inverted  for  the  model  coefficients  of  global  phase  velocity 
maps  expanded  in  spherical  harmonics.  The  predicted  phase  anomaly  between  a  source  at 


4 


(^5)  (t>s)  and  a  receiver  at  {Or,  is  written  in  terms  of  real  spherical  harmonics 

= - /  Y^Y{Ai^cosm<i)  +  Bim5mm(t>)p'P{e)ds,  (1) 

Co  J(es,(l>s)  1=0  m=0 

where  Aim,  Bim  are  the  model  coefficients,  L  is  the  maximum  angular  degree  of  the  ex¬ 

pansion,  and  are  the  Legendre  polynomials.  The  coefficients  Aim,  Bim  are  determined 
by  damped  least-squares  inversion.  The  result  of  these  inversions  are  global  phase  velocity 
maps  {Ekstrom  et  al,  1997).  The  model  coefficients  can  be  obtained  in  electronic  form  via 
anonymous  ftp  to  saf.harvard.edu,  by  retrieving  the  files  in  pub/ETL-JGR-96. 

The  resulting  maps  explain  80-96%  of  the  variance  in  Love  wave  phase  anomalies,  and 
between  70-90%  of  the  Rayleigh  wave  data.  The  lowest  variance  reduction  is  obtained  for 
the  longest  periods.  Figure  1  shows  the  phase  velocity  map  for  Eurasia  and  Love  waves  at 
35  s  period.  Clearly,  a  very  large  portion  of  the  signal  is  associated  (as  expected)  with  the 
contrast  between  oceanic  and  continental  crust  and  lithosphere.  Within  Eurasia,  very  large 
variations  are  seen  which  reflect  both  the  total  thickness  of  the  crust  (as  beneath  western 
Tibet),  tectonic  activity  (the  Alpine-Himalayan  belt  in  general)  and  very  thick  sedimentary 
basins  (.\rctic  ocean  shelf,  southern  Caspian  sea).  Resolution  experiments  indicate  that  these 
features  are  well  resolved.  It  is  interesting  to  note  that  when  compared  with  phase  velocity 
predictions  from  the  recent  global  crustal  model  CRUST-5.0  {Mooney  et  al.,  1995),  there  are 
some  areas  which  show  very  large  similarities,  and  other  areas  (in  particular  sedimentary 
basins)  which  do  not.  The  correlation  between  the  observed  phase  velocity  map  and  that 
predicted  by  Mooney  et  al.'s  model  is  approximately  0.80  for  Love  waves  at  35  s. 

A  degree  20  model  of  the  Mantle 

The  surface  wave  data  described  above  have  been  u.sed  in  conjunction  with  previous  seis- 
mological  datasets  to  derive  a  higher  resolution  3-D  model  of  the  Earth’s  mantle.  Prelimi¬ 
nary  results  have  been  presented  at  various  nuH'tings  {Ekstrom  and  Dziewonski,  1994;  1995; 
Dziewonski  et  al.,  1996). 

The  parameterization  and  inversion  approach  are  similar  to  those  used  for  earlier  Harvard 
mantle  models,  such  as  S12WM13  (Su  et  al..  1993).  with  the  following  notable  differences: 

•  We  use  a  split  parameterization  for  the  mantle,  aad  parameterize  the  radial  variations 
in  the  upper  mantle  in  terms  of  Chebyshev  polynomials  of  degree  0-7.  For  the  lower 
mantle  we  use  Chebyshev  polynomials  of  degree  0-5. 

•  The  spherical  harmonic  expansion  is  exteiuhnl  from  degree  12  to  degree  20. 

•  The  data  are  corrected  for  crustal  structure  using  a  new  compilation  of  crustal  thick¬ 
nesses  and  velocities  by  Mooney  et  al.  (1995).  This  is  in  contrast  to  previous  models 
which  were  only  corrected  for  average  continental  and  oceanic  crustal  structure. 

Several  types  of  data  with  sensitivity  to  velocity  variations  in  different  parts  of  the  mantle 
are  included.  The  following  data  have  been  incorporated  in  the  derivation  of  the  new  model 

•  Absolute  and  differential  travel  times  measured  from  waveform  data  {S,  SS,  ScS, 
SS  -  5,  ScS  -  S,  ScS  ScS  -  ScS.  S  -  SKS,  SKKS  -  SKS). 

•  Long-period  body  and  mantle  waveforms  from  earthquakes  recorded  on  the  Global 
Seismic  Network. 
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Figure  2:  Map.s  .show  the  5  velocity  variations  at  100  km  depth  beneath  the  Pacific  and 
adjacent  areiis.  The  top  panel  shows  the  S-velocity  variations  recovered  from  primarily 
P/S\'-sensitive  data,  and  the  bottom  panel  for  SH-sen.sitive  data.  There  is  strong  similarity 
between  the  two  maps  in  most  regions. 


•  \’ery  long-period  mantle  waveforms  from  large  earthquakes  recorded  on  the  Global 
Seismic  Network. 

•  Intermediate  period  dispersion  measurements  of  Rayleigh  and  Love  waves  (described 
above). 

The  new  model  S20U7L5  shows  many  of  the  same  features  seen  in  previous  tomographic 
images,  but  with  significantly  sharper  definition.  On  a  global  scale,  the  fast  velocities  beneath 
stable  continental  interiors  are  the  most  prominent  anomalies  in  the  top  300  km  of  the  mantle. 
The  very  high  velocities  beneath  the  stable  interior  of  West  Africa  continue  even  deeper.  The 
new  model  also  shows  that  the  slow  velocity  anomalies  associated  with  ridges  continue  to 
great  depth,  in  particular  in  the  Indian  Ocean.  The  Red  Sea  rift  shows  slow  velocities 
extending  down  to  400  km.  Below  250  km  slow  features  not  directly  associated  with  surface 
tectonics  predominate,  in  particular  in  the  center  of  the  Pacific  Ocean. 

In  order  to  investigate  the  compatibility  of  the  various  datasets,  we  have  performed  several 
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Figure  3:  Maps  show  the  S  velocity  variations  at  200  km  depth  beneath  the  Pacific  and 
adjacent  areas  obtained  from  separate  inversions  of  P/SV  and  SH  sensitive  data.  Note  the 
slow  velocities  in  the  central  Pacific  seen  in  the  top  (P/SV)  map,  but  absent  in  the  bottom 
(SH)  map.  .Also,  the  East  Pacific  Rise  remains  a  slow  feature  at  this  depth  for  SH  data,  but 
not  for  P/S\’  data. 


inversions  using  subsets  of  the  full  data.set.  In  particular,  we  have  separated  the  data  which 
are  primarily  sensitive  to  SH-type  structure  from  those  which  primarily  sense  P/SV-type 
structure,  by  considering  separately  observations  made  from  transversely  polarized  seismo¬ 
grams  and  those  made  on  vertical  and  radial  sei.smograms.  The  somewhat  surprising  result 
is  that  there  are  regions  and  depth  ranges  where  we  obtain  distinctly  different  patterns  in 
the  separate  inversions.  In  Figure  2  we  show  a  comparison  of  the  structures  at  100  km  depth 
beneath  the  Pacific  as  imaged  using  the  two  sets  of  data.  While  there  are  differences  between 
the  two  maps,  overall  there  is  fairly  close  agreement.  Mostly  we  see  the  fast  roots  beneath 
the  continents,  and  the  signature  of  the  hot  and  cooling  lithosphere  in  the  oceans.  At  200  km 
(Figure  3),  however,  the  slow  velocity  in  the  central  Pacific  seen  in  the  P/SV  map  (as  well 
as  in  many  previous  models)  is  not  seen  in  the  SH  map,  which  instead  shows  consistently 
fast  velocities  in  the  same  region.  The  SH  map  also  shows  slower  velocities  at  this  depth 
associated  with  the  East  Pacific  Rise,  while  the  P/SV  map  does  not. 
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There  are  different  possible  explanations  for  this  result.  The  most  straightforward  reason 
could  be  that  this  is  the  effects  of  varying  levels  of  transverse  anisotropy  in  the  asthenosphere. 
A  second  possibility  is  that  azimuthal  anisotropy  in  the  oceanic  lithosphere  and/ or  astheno¬ 
sphere  gets  mapped,  through  the  ray  coverage,  into  different  apparent  velocity  structures  for 
Love  and  Rayleigh  waves.  A  third  possibility  is  that  we  have  overestimated  the  resolving 
power  of  the  different  datasets,  and  that  the  observed  patterns  are  spurious.  Based  on  the 
robustness  of  the  results  (established  through  additional  experiments),  we  do  not  believe 
the  third  explanation  is  likely,  and  we  are  currently  exploring  the  need  for  including  specific 
considerations  of  anisotropy  in  future  modeling  efforts. 

Modeling  of  Regional  Phases 

The  limited  distribution  of  sources  and  the  sparse  (in  a  global  sense)  distribution  of  re¬ 
ceivers  set  a  limit  on  our  ability  to  simply  keep  increasing  the  resolution  of  our  global  models. 
To  obtain  good  global  resolution  on  the  scale  of  500  km  or  less  will  probably  be  very  difR- 
cult,  if  not  impossible,  for  some  time  to  come.  It  is  possible  though  to  obtain  high  resolution 
models  of  a  number  of  regions,  both  large  and  small,  that  are  well  sampled.  To  obtain  the 
best  possible  coverage  in  a  given  region,  one  would  ideally  use  sources  and  receivers  that  are 
both  within  and  outside  the  region  of  interest.  Proper  treatment  of  such  a  combination  of 
data  would  require  a  hybrid  parameterization  {Fukao  et  al,  1992). 
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Figure  4:  Map  showing  ray  coverage  for  Pn  observations  in  northeastern  U.S.,  reported  in 
the  ISC. 

In  a  preliminary  investigation  we  have  been  examining  how  well  our  current  3-D  mantle 
models  may  be  used  to  explain  trends  in  regional  data  in  different  regions.  In  Figure  4  we 
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show  the  Pn  raypath  coverage  for  events  reported  by  the  ISC  for  northeastern  North  Aiherica. 
This  region  was  chosen  to  test  our  model  in  an  area  where  we  predict  large  anomalies  based 
on  the  model  S&:P12WM13  {Su  and  Dziewonski,  1993)  and  have  a  good  path  coverage  using 
data  from  the  ISC. 


Residuals  Predicted  using  S&P12/WM13  (s) 

Figure  5:  Observed  residuals  for  Pn  compared  with  predictions  from  S&;P12WM13  for  paths 
in  northeastern  U.S. 

Figure  5  shows  the  observed  and  calculated  Pn  residuals  for  this  region.  The  Pn  travel 
time  was  calculated  using  the  average  P  velocity  at  the  top  of  the  mantle  in  the  model 
Si:P12\VM13.  Despite  some  scatter  a  correlation  between  the  residuals  found  for  these 
regional  pha,scs  and  those  derived  from  our  global  model  can  be  seen.  In  particular,  the 
level  of  heterogeneity  appears  to  be  well  matched.  A  second  region  for  which  we  performed 
this  experiment  was  southern  Eurasia.  This  region  has  many  fewer  stations  reporting  to  the 
ISC.  Figure  6  shows  the  ISC  residuals  and  the  predicted  residuals.  In  this  case  there  is  no 
apparent  correlation.  This  may  be  for  one  of  several  reasons:  there  is  some  evidence  that 
earthquakes  in  this  region  are  poorly  located:  our  model  in  this  region  may  not  be  so  well 
constrained  due  to  lack  of  data  coverage:  structure  in  this  region  may  be  more  complex  than 
the  resolution  of  our  current  model  allows;  or  the  correlation  may  be  obscured  due  to  the 
smaller  range  of  predicted  anomalies. 

Despite  being  based  wholly  on  teleseismic  data,  our  current  global  models  appear  to  pro¬ 
vide  fairly  good  estimates  of  uppermost  mantle  velocities  beneath  regions  such  as  northeast¬ 
ern  North  America.  The  experiment  also  suggests  that  beneath  areas  such  as  Eurasia  regional 
phases  may  need  to  be  included  to  adequately  model  the  uppermost  mantle  structure.  How¬ 
ever,  such  an  analysis  will  require  careful  relocation  of  the  seismicity  before  utilizing  the 
travel  times.  With  the  inclusion  of  regional  data,  parameterization  in  terms  of  spherical 
harmonics  is  a  poor  choice  since  by  definition  they  have  a  uniform  resolution  everywhere.  In 
future  work,  we  will  instead  move  towards  hybrid  parameterizations  with  variable  resolution, 
to  be  able  to  include  and  subsequently  to  better  model  observations  corresponding  to  both 
regional  and  teleseismic  distances. 
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Residuals  Predicted  using  S&P12/WM13  (s) 

Figure  C;  Observed  residuals  for  Pn  compared  with  predictions  from  S&;P12WM13  for  paths 
in  southern  Eurasia. 


Conclusions  and  Recommendations 

By  combining  a  large  number  of  diverse  seismological  observations,  we  have  made  progress 
towards  obtaining  higher  resolution  tomographic  images  of  the  regional  scale  elastic  structure 
globally  and  beneath  Eurasia.  The  shallowest  part  of  the  mantle  is  effectively  constrained  and 
imag('(l  using  intermediate  period  surface  waves.  Both  global  and  regional  scale  heterogeneity 
in  the  mantle  has  been  shown  to  bias  seismically  derived  event  locations.  In  some  regions, 
our  current  tomographic  models  can  predict  the  gross  residual  patterns  seen  in  regional  phase 
data.  It  appears  likely  that  further  improvements  in  locations  can  be  achieved  by  correcting 
for  better  known  elastic  heterogeneity  in  the  Earth's  mantle. 
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Abstract 

A  new  technique  for  making  single-station  phase  velocity  measurements  is  developed  and  applied 
to  a  large  number  of  globally  recorded  Rayleigh  and  Love  waves  in  the  period  range  35-150  s. 
The  method  is  based  on  phase-matched  filter  theory  and  iteratively  suppresses  the  effect  of 
interfering  overtones  by  minimizing  residual  dispersion.  The  model  surface  wave  signal  is 
described  by  its  amplitude  and  apparent  phase  velocity,  both  of  which  are  parameterized  in 
terms  of  smooth  B-spline  functions  of  frequency.  A  misfit  function  is  constructed  which 
represents  the  difference  between  the  model  and  observed  waveforms,  and  the  optimal  spline 
coefficients  are  estimated  in  an  iterative  misfit  minimization  algorithm.  In  order  to  eliminate 
cycle  skips  in  the  measurements  of  phase  at  short  periods,  the  waveforms  are  first  matched  at 
long  periods,  and  the  frequency  range  is  gradually  extended  to  include  higher  frequencies.  The 
application  of  the  algorithm  to  records  from  the  Global  Seismographic  Network,  using 
earthquakes  in  the  Harvard  centroid-moment  ten.sor  catalog,  results  in  the  determination  of  more 
than  50.000  high-quality  dispersion  curves.  The  observed  variations  in  measured  dispersion  for 
pairwise  similar  paths  are  used  to  estimate  realistic  uncertainties  in  the  data.  Phase  delays  at 
discrete  periods  are  inverted  for  global  maps  of  variations  in  phase  velocity  expanded  in  spherical 
harmonics  up  to  degree  40.  A  realistic  resolution  test  indicates  that  structures  are  well  recovered 
up  to  at  least  degree  20.  The  new  phase  velocity  maps  explain  70-96%  of  the  observed  variance 
in  phase  residuals,  reflecting  the  high  internal  consistency  of  the  dispersion  measurements. 
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Introduction 

Regional  variations  in  the  propagation  speeds  of 
teleseismic  surface  waves  have  been  noted  since  the 
early  part  of  the  century  [Tams,  1921a,  b].  In  the 
1960s,  the  magnitudes  of  these  global  variations  for 
weaves  in  the  period  range  1O-150  s  could  be  summa¬ 
rized  and  discussed  in  terms  of  gross  lateral  differences 
in  the  layering  and  structure  of  the  crust  and  up¬ 
per  mantle  [e.g.,  Oliver,  1962;  Dorman,  1969].  Vari¬ 
ations  at  longer  periods  were  also  observed  and  in¬ 
terpreted  in  terms  of  lateral  variations  of  upper  man¬ 
tle  elastic  and  anelastic  structure  [e.g.,  Toksdz  and 
Anderson,  1966;  Kanamori,  1970;  Dziewonski,  1971; 
Wu,  1972].  With  the  start  of  digital  data  collec¬ 
tion  from  new  global  seismograph  networks  [Agnew 
et  al,  1976;  Peterson  et  al,  1976]  and  the  develop¬ 
ment  of  new  techniques  for  analyzing  surface  wave 
dispersion  through  waveform  inversion  [e.g.,  Dziewon¬ 
ski  and  Steim,  1982],  efforts  began  to  map  the  global 
variations  in  long-period  phase  velocities  and  the  cor¬ 
responding  three-dimensional  (3-D)  structure  of  the 
Earth  using  tomographic  techniques  [Masters  et  al, 
1982;  Nakanishi  and  Anderson,  1982,  1983,  1984; 
Woodhouse  and  Dziewonski,  1984;  Tanimoto  and  An¬ 
derson,  1985;  Nataf  et  al.,  1986;  Tanimoto,  1986; 
Wong,  1989;  Montagner  and  Tanimoto,  1990,  1991]. 
These  studies  have  made  use  of  a  variety  of  seismic 
data  to  image  the  Earth’s  interior;  surface  waves  have, 
however,  provided  the  most  valuable  constraints  on 
upper  mantle  structure. 

Although  good  recordings  of  intermediate-period 
surface  waves  are  much  more  common  than  good 
very  long-period  recordings,  the  shortest-period  sur¬ 
face  waves  included  in  the  global  tomographic  studies 
referenced  above  is  75  s,  and  most  of  the  studies  only 
considered  surface  waves  with  periods  greater  than 
150  s.  Two  explanations  for  this  situation  can  be 
found.  First,  while  shorter-period  surface  waves  are 
sensitive  to  velocity  variations  in  the  shallow  mantle, 
they  also  depend  critically  on  the  structure  and  thick¬ 
ness  of  the  crust,  which  are  not  very  w*ell  known  on  the 
global  scale.  By  limiting  surface  wave  observations  to 
periods  which  have  less  sensitivity  to  the  crust  and 
by  applying  simple  a  priori  crustal  corrections  to  the 
data,  these  studies  have  sought  to  eliminate  crustal 
effects.  Second,  it  is  more  difficult  to  measure  the 
dispersion  of  surface  waves  at  intermediate  periods 
than  at  long  periods.  In  particular,  establishment  of 
the  total  propagation  phase  of  a  teleseismic  surface 
wave  at  a  particular  period,  such  that  there  is  no  re¬ 
maining  ambiguity  in  the  integral  number  of  cycles, 


is  a  nontrivial  task. 

There  axe,  in  addition,  several  complicating  as¬ 
pects  of  surface  wave  propagation  in  a  laterally  het¬ 
erogeneous  Earth  which  are  more  significant  for  shorter- 
period  waves,  but  which  are  not  very  well  understood 
or  predicted,  such  as  focusing,  refraction  of  the  ap¬ 
parent  ray  path,  and  scattering  [e.g.,  Evemden,  1954; 
Capon,  1970;  Sobel  and  von  Seggem,  1978;  Lay  and 
Kanamori,  1985;  Woodhouse  and  Wong,  1986',  Laske, 
1995].  These  complications  usually  do  not  prevent  us 
from  making  measurements  of  dispersion  but  do  sug¬ 
gest  that  the  interpretation  of  these  measurements  in 
terms  of  Earth  structure  will  depend  on  whether  the 
chosen  theoretical  model  for  surface  wave  propagation 
in  a  heterogeneous  Earth  is  adequate. 

Recently,  interest  in  incorporating  shorter-period 
surface  waves  in  global  studies  has  grown  [Trampert 
and  Woodhouse,  1995,  1996;  Zhang  and  Lay,  1996; 
Laske  and  Masters,  1996;  Ekstrom  and  Dziewonski, 
1995].  This  is  primarily  due  to  the  high  resolving 
power,  both  radially  and  laterally,  that  these  waves 
can  afford  in  the  shadowiest  mantle.  Careful  investi¬ 
gation  of  these  waves  can  help  resolve  some  difficult 
geodynamical  questions,  such  as  the  issue  of  “pas¬ 
sive”  (shallow)  versus  “active”  (deep)  mechanisms  for 
spreading  at  mid-ocean  ridges  [Zhang  and  Tanimoto, 
1992;  Su  et  al,  1992]  and  the  depth  extent  and  chem¬ 
ical  versus  thermal  nature  of  lithospheric  roots  under 
parts  of  the  continents  [Jordan,  1975]. 

In  addition  to  providing  information  about  the 
elastic  and  anelastic  structure  of  the  crust  and  shallow 
mantle,  intermediate-period  surface  waves  are  impor¬ 
tant  for  earthquake  source  studies.  For  shallow  earth¬ 
quakes,  surface  waves  are  commonly  the  largest  tele¬ 
seismic  phases  recorded  on  long-period  seismograms. 
For  recordings  with  a  low  signal- to-noise  ratio,  they 
may  be  the  only  part  of  the  seismogram  that  can 
be  used  for  a  source  mechanism  determination.  Im¬ 
provements  in  our  ability  to  model  global  surface  wave 
propagation  at  shorter  periods  will  make  it  possible 
to  routinely  derive  source  geometries  and  moments  for 
earthquakes  in  the  magnitude  range  4.5  <  M  <  5.0, 
something  which  is  now  feasible  only  in  areas  where 
there  is  access  to  regional  or  local  data.  It  is  worth 
noting  that  for  this  application,  the  interpretation  of 
the  dispersion  in  terms  of  a  3-D  Earth  model  is  of  less 
importance. 

In  this  paper,  our  first  objective  is  to  present  a 
new  analysis  algorithm  which  overcomes  the  difficulty 
of  determining  the  absolute  propagation  phase  for 
highly  dispersed  intermediate-period  surface  waves. 
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Our  second  objective  is  to  present  the  data  set  of  dis¬ 
persion  measurements  and  associated  realistic  uncer¬ 
tainties  which  has  resulted  from  the  systematic  appli¬ 
cation  of  the  new  algorithm  to  a  very  large  volume 
of  global  digital  data.  Our  final  objective  is  to  derive 
and  present  phase  velocity  maps,  which  can  be  used 
to  model  global  surface  wave  propagation  in  the  pe¬ 
riod  range  35-150  s,  as  well  as,  in  future  studies,  to 
infer  the  3-D  structure  of  the  shallow  mantle. 

Theory 

We  describe  and  interpret  the  propagation  of  fun¬ 
damental  Love  and  Rayleigh  waves  using  ray  theory 
on  a  sphere  [e.g.,-  Tromp  and  Dahlen^  1992,  1993]  and 
write  the  surface  wave  seismogram  u{uj)  as 

=  A{uj)  exp[i^(a;)],  (1) 

where  and  are  the  amplitude  and  phase, 
respectively,  of  the  wave  as  functions  of  the  angular 
frequency  u;.  For  a  given  source-receiver  geometry, 
the  phase  ^  is  the  sum  of  four  terms, 

^  =  $5  +  -h  (2) 

where  ^5  is  the  source  phase  calculated  from  the 
source  mechanism  and  geometrical  ray  takeoff  az¬ 
imuth,  is  the  receiver  phase,  is  the  static  phase 
contribution  from  each  ray  focus,  and  is  the  prop¬ 
agation  phase 

where  c  is  the  local  phase  velocity  and  the  integra¬ 
tion  follows  the  ray  path.  The  amplitude  -4  can  be 
expressed  as 

.4  =  .45-4p.4j^.4q,  (4) 

where  .45  is  the  excitation  at  the  source,  Ar  is  the 
receiver  amplitude,  .4a  is  the  geometrical  spreading 
factor,  and  Aq  is  the  decay  factor  due  to  attenua¬ 
tion  along  the  ray  path.  When  the  location  and  focal 
mechanism  of  the  earthquake  are  known,  a  theoreti¬ 
cal  reference  seismogram  u^(u;)  based  on  a  spherical 
Earth  model  can  be  calculated  and  written  as 


where  c®  is  the  spherical  Earth  phase  velocity,  A  is 
the  angular  epicentraJ  distance,  R  is  the  radius  of  the 
Earth,  and  X  is  the  propagation  path  length  mea¬ 
sured  along  the  great  circle.  We  express  the  observed 
surface  wave  u{u)  as  a  perturbation  with  respect  to 
the  reference  seismogram 

u(u/)  =  [A°(a;)  +  M(w)]  expi[$°(a;)  +  (54>(a;)].  (7) 

Attributing  to  a  perturbation  of  the  propagation 
phase,  we  have 

=  +  =  (8) 

-h  dc 

where  5c  is  the  apparent  average  phase  velocity  per¬ 
turbation,  calculated  for  the  distance  X  along  the 
great  circle. 

In  order  to  find  an  adequate  and  simple  model  pa¬ 
rameterization  for  a  realistic  observed  surface  wave 
shape,  we  assume  that  the  phase  and  amplitude  of 
Love  and  Rayleigh  waves  vary  smoothly  with  fre¬ 
quency  and  that  discontinuities  in  the  phase,  which 
can  result,  for  example,  from  multipathing,  are  rare. 
Evidence  for  refraction  and  multipathing  at  periods 
shorter  than  30  s  is  not  uncommon  [e.g.,  Lemer-Lam 
and  Park^  1989;  Levshin  et  aL,  1992],  but  our  experi¬ 
ence  is  that  a  smooth  parameterization  of  phase  and 
amplitude  is  sufficient  to  model  most  of  the  observa¬ 
tions.  We  therefore  express  the  amplitudes  as 

N 

(9) 

where  6,(z  =  1,2, ...,  A^)  are  spline  coefficients  for  the 
cubic  B-spline  polynomials  fi  (see  Figure  1)  and  bo  is 
a  constant.  Initially,  the  coefficients  bi  are  calculated 
for  a  reference  Earth  model  and  normalized  such  that 

=  (10) 

1=1 

The  propagation  phase  is  a  rapidly  growing 
function  of  lj  and  therefore  not  eatsily  splined.  In¬ 
stead,  we  spline  the  apparent  average  phase  velocity 
perturbation  Sc 


u^{u;)  =  .4°(a;)exp[2^^(Lt;)]. 


(5) 


The  propagation  phase  for  the  reference  surface  wave 
is 


ijRA 


u;A 


(6) 


N 

^c{u>)  =  'Y^difi{u)),  (11) 

1=1 

where  fi  are  the  same  cubic  B-spline  polynomials  as 
above  and  di  are  spline  coefficients  to  be  determined 
for  each  path. 
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This  parameterization  will  not  be  adequate  for  all 
observations.  For  example,  for  paths  along  which 
multipathing  or  scattering  is  severe,  we  will  not  achieve 
a  good  agreement  between  the  observed  seismogram 
and  this  model  of  the  signal.  However,  by  applying  a 
selection  criterion  bcised  on  fit  between  the  observed 
and  model  signal,  we  can  identify  and  eliminate  many 
of  these  paths  and  corresponding  unreliable  dispersion 
measurements. 

Measurement  Technique 

The  objective  is  to  separate  out  the  fundamen¬ 
tal  mode  surface  wave  in  the  seismogram  and  to 
make  a  measurement  of  its  dispersion.  Real  data 
contain  noise,  as  well  as  interfering  seismic  phases, 
so  we  nec‘d  a  method  which  suppresses  noise  and 
isolat(»s  the  fundamental  mode  surface  waves  from 
oth(‘r  phases.  Since  we  will  require  many  measure¬ 
ments.  w(*  also  need  a  method  which  can  be  auto- 
mat(^d.  Many  techniques  have  been  presented  for 
making  singh^station  dispersion  measurements;  the 
method  we  develop  builds  on  the  ideas  of  residual  , 
disj)ersion  [Dztrwonski  and  Hales,  1972;  Dziewonski 
et  ai,  1972]  and  phase-matched  filtering  [e.g.,  Herrin 
and  Goforth.  1977).  Our  method  adds  several  features 
to  the.s(*  two  basic  concepts. 

First,  we  cxssume  that  the  observed  seismogram 
contains  a  dispersed  surface  wave  signal  .4^(u;)  exp[/4>‘'(-^’)]‘ 
wh(T(*  the  S  superscript  refers  to  the  observed  sig¬ 
nal.  The  phase'- matched  filter  which  will  optimally 
enhance  and  compress  this  signal  in  the  presence  of 
white  noise*  is  the  signal  it.self  (e.g.,  Herrin  and  Go¬ 
forth.  1977].  That  is.  by  autocorrelating  the  signal, 
the  dispersion  is  annihilated  and  the  noi.se  is  su{>- 
pressed:  the  result  is  a  zero-phase  filtered  spike  with 
amplitude  spectrum 

Second,  we  consider  how  to  best  isolate  the*  funda¬ 
mental  mode*  from  interfering  overtones.  This  is  also 
achieved  by  correlation  with  the  phase-annihilating 
filter  exp[?4>'''(,i^’)].  However,  only  the  part  of  the* 
time  domain  correlation  function  close  to  its  nuixi- 
mum  corresponds  to  fundamental  mode  energy.  In 
order  to  isolate  the  fundamental  mode  energy.  th(* 
part  of  the  correlation  function  near  the  imiximum 
can  be  windowed  in  time  [e.g.,  Dziewonski  et  at.,  1972; 
Herrin  and  Goforth,  1977;  Lemer-Lam  and  Jordan, 

1983;  Levshin  et  al,  1992].  If  the  phase  delay  is  then 
reintroduced  by  multiplying  the  windowed  correlation 
function  by  exp[4>^(j;)]  in  the  spectral  domain,  a  fil¬ 
tered  version  of  the  input  signal  is  obtained,  in  which 


signals  with  dispersion  different  from  that  of  the  fun¬ 
damental  mode  are  suppressed.  However,  there  is  no 
single  time  window  which  will  be  optimal  for  win¬ 
dowing  the  correlation  function.  If  a  broad  range  of 
frequencies  is  considered  and  the  window  is  selected 
to  include  a  cycle  or  more  of  the  longest  periods,  then 
the  window  used  will  be  too  wide  for  high  frequencies 
(thereby  including  signals  that  have  different  group 
delays)  and  vice  versa. 

On  the  basis  of  these  considerations,  we  develop 
the  following  algorithm  for  estimating  the  dispersion 
and  amplitude  of  the  fundamental  mode  surface  wave 
in  the  observed  signal.  A  trial  fundamental  mode 
model  seismogram  A^{u)  exp[i(t>^^  {uj)]  is  constructed 
based  on  the  source  and  receiver  locations,  the  focal 
mechanism,  and  the  phase  and  amplitude  effects  of 
propagation  in  the  reference  Earth  model.  We  then 
construct  the  whitening  phase-matched  filter  W 

W{u) 

and  cross  correlate  this  function  with  the  observed 
seismogram.  If  the  observed  seismogram  contained  a 
signal  with  identical  dispersion  and  amplitude  spec¬ 
trum  to  our  trial  synthetic  seismogram  and  if  we  were 
considering  all  frequencies,  the  result  of  the  cross  cor¬ 
relation  would  be  a  delta  function  in  time.  If  we  con¬ 
sidered  a  range  of  frequencies,  the  result  would  be  a 
band-passed  delta  function.  The  trial  dispersion  and 
amplitude  will,  in  general,  not  match  the  observed 
signal,  and  there  will  be  a  misfit.  By  adjusting  the 
model  seismogram  to  minimize  the  misfit,  an  opti¬ 
mal  estimate  of  the  dispersion  and  amplitude  in  the 
observed  signal  is  obtained. 

The  specific  misfit  function  must  be  chosen  with 
care,  and  we  have  found  the  following  approach  to 
b(‘  efficient  and  useful.  Consider  the  full  range  of 
frc*(iuencies  [flmin,  ^max]  for  which  a  dispersion  and 
amplitude  estimate  is  to  be  obtained.  In  this  study, 
the  full  period  range  is  250-32  s.  Divide  this  range 
into  2iV  —  1  subranges  such  that 

^min  ^  f^min  +  “  l)Aa;,  .  . 

(^rnax  —  ^min  +  (^  “h  l)Au;, 

for  i  —  1,2, ...,  2N  —  1  and  Au;  =  (f^max  ““ Dmin)/ (2A^), 
and  construct  2N  —  1  narrow  band-pass  filters  Fj(a;) 
using  Welch  tapers  [Press  et  a/.,  1986]  over  these 
ranges  (Figure  2).  Filter  the  cross-correlation  func¬ 
tion  in  each  of  these  bands, 


Giioj)  Fi{u)S{uj)W^{uj) 


(14) 
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and  do  the  same  to  the  whitened  autocorrelation  of 
the  trial  synthetic  seismogram, 

Hiiuj)  =  Fi{Lj)A^{uj)exp[i^^{uj)]W'{Lj) 

=  Fiiu)  ,  (15) 


that  is,  a  band-passed  delta  function.  These  functions 
are  Fourier  transformed  back  to  the  time  domain  to 
obtain  the  band-passed  cross  correlation  gi{t)  and  the 
model  autocorrelation  hi{t)  to  form  the  misfit  func¬ 
tion 


2A^-l 


Wi 


t=l 


ZjL-jMtjW 


(16) 


where  Ck  and  d^.  are  the  spline  coefficients  to  be  v'aried, 
iL't  is  a  weight  given  to  the  misfit  in  the  zth  frequency 
passband,  and  Ji  is  the  width  (in  terms  of  time  point 
samples)  of  the  time  window*  over  which  the  cross  cor¬ 
relations  and  autocorrelations  are  compared.  We  gen¬ 
erally  choose  Ji  to  include  three  cycles  of  the  center 
period  in  the  passband,  that  is,  Ji  =  3T/(2Af),  wdiere 
T  is  the  center  period  47r/(u;J„^x 
sampling  interval. 

The  main  difficulty  in  obtaining  dispersion  mea¬ 
surements  from  teleseismically  observed  surface  w’aves 
at  periods  less  than  about  100  s  is  that  variations  in 
phase  velocity  are  sufficiently  large  to  cause  path  vari¬ 
ations  of  several  cycles  in  phase,  and  isolated  residual 
phase  measurements  in  the  range  [-7r,7r]  at  a  single 
period  are  therefore  ambiguous  in  multiples  of  2;r.  For 
minor  arc  observations,  how'ever,  phcise  measurements 
at  periods  longer  than  100  s  can  be  associated  with 
a  total  cycle  count  without  difficulty.  If  a  continuous 
dispersion  curve  can  be  anchored  at  long  periods,  the 
total  phase  perturbation  can  then  be  inferred  without 
ambiguity. 

In  order  to  obtain  the  complete  broadband  disper¬ 
sion.  we  therefore  employ  a  method  w'e  call  iterative 
frequency  band  expansion  in  the  minimization  of 
We  initially  determine  the  dispersion  parameters  in  a 
narrow  frequency  range  centered  around  100  s  and 
then,  in  several  iterations,  increase  the  higher  fre 
quency  end  of  the  total  frequency  range.  In  each 
iteration,  the  dispersion  and  amplitude  parameters 
are  adjusted  to  minimize  ^  using  the  downhill  sim¬ 
plex  method  of  Nelder  and  Mead  [1965].  Only  the 
frequency  bands  which  fall  within  the  range  of  con¬ 
sidered  frequencies  are  included  in  the  misfit  function 
sum  (equation  (16)).  A  search  algorithm  such  as  the 
dowmhill  simplex  method  is  convenient  for  solving  this 
problem,  since  it  is  easy  to  incorporate  constraints 


such  as  nonnegativity  of  A  and  realistic  limits  on  Sc, 
Once  ^  is  minimized  for  a  given  frequency  range,  the 
range  is  expanded  to  include  higher  frequencies.  Be¬ 
cause  of  the  smoothness  of  Sc,  a  prediction  can  be 
made  for  frequencies  slightly  higher  than  the  ones  in¬ 
cluded  in  the  previous  iteration.  This  ensures  that 
the  phase  difference  between  model  seismogram  and 
observed  seismogram  always  remains  small,  so  that 
gradient  methods  for  minimizing  ^  continue  to  w*ork 
in  the  subsequent  iteration. 

An  example  of  the  minimization  procedure  is  given 
in  Figure  3.  The  observed  seismogram,  filtered  be¬ 
tween  250  and  32  s,  is  shown  together  with  a  model 
seismogram  calculated  for  the  combination  of  the 
mantle  heterogeneity  model  SH8U4L8  [Dziewonski 
and  Woodward,  1992]  and  the  preliminary  reference 
Earth  model  (PREM)  [Dziewonski  and  Anderson, 
1981].  The  two  waveforms  are  quite  different,  and  it 
is  clear  that  major  adjustments  are  required  in  both 
phase  and  amplitude  in  order  to  make  them  agree. 
When  filtered  between  250  and  75  s  (Figure  3b)  there 
is  better  agreement.  After  adjustment  of  phase  and 
amplitude  in  this  period  range,  two  w^aveforms 
are  very  similar  (Figure  3c);  the  frequency  window  is 
broadened  (Figure  3d);  again,  the  phase  and  ampli¬ 
tude  are  adjusted.  After  five  iterations  of  frequency 
w'indow’  expansion,  the  full  waveform  is  matched  (Fig¬ 
ure  3e).  The  resulting  dispersion  and  amplitude  per¬ 
turbations  are  showm  in  Figure  4.  No  damping  of  the 
dispersion  curve  to  the  starting  curve  is  involved,  and 
the  results  are  largely  insensitive  to  the  initial  values 
of  .4  and  Sc. 

Application 

We  have  applied  the  method  described  in  the  pre¬ 
vious  section  to  digital  seismograms  from  the  Global 
Seisinographic  Network  (GSN)  of  Incorporated  Re¬ 
search  Institutions  for  Seismology  (IRIS),  the  Chinese 
Digital  Seismograph  Netw^ork  (CDSN),  the  Global 
Telemetered  Seismograph  Network  (GTSN),  and  the 
MEDNET  and  GEOSCOPE  networks,  using  earth¬ 
quake  sources  from  the  period  January  1989  to  Septem¬ 
ber  1995.  Figure  5  shows  the  distribution  of  the  earth¬ 
quakes  and  stations  from  which  data  have  been  col¬ 
lected. 

Analysis  w*as  attempted  for  earthquakes  with  Mw  > 
5.5.  We  used  the  moment  tensors  and  centroid  loca¬ 
tions  in  the  Harvard  centroid-moment  tensor  (CMT) 
catalog  [Dziewonski  et  aL,  1981]  to  calculate  the 
source  excitation.  Only  shallow  earthquakes  {h  < 
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50  km)  were  included  in  order  to  maximize  the  ex¬ 
citation  of  surface  waves.  In  addition,  we  considered 
only  paths  in  the  distance  range  25®  <  A  <  150°  in 
order  to  avoid  problems  inherent  in  isolating  the  fun¬ 
damental  mode  at  short  distances  and  close  to  the 
antipode.  The  variations  in  5c  and  A  between  250  s 
and  32  s  were  parameterized  using  six  spline  basis 
functions.  After  several  experiments,  this  number  of 
spline  bmctions  was  found  to  be  adequate  for  describ¬ 
ing  most  of  the  variations  seen  in  the  data. 

For  each  seismogram,  the  instrument  response  was 
removed  to  restore  ground  displacement.  Rayleigh 
wave  dispersion  was  determined  from  the  vertical 
component,  and  Love  wave  dispersion  was  determined 
from  transverse  records  constructed  by  rotation  of  the 
horizontal  components  using  the  great  circle  back  az¬ 
imuth.  More  than  128,000  station-receiver  paths  were 
analyzed,  but  only  a  fraction  of  these  yielded  useful 
meiisiirements.  For  some  paths,  we  obtained  good 
dispersion  curves  to  45-50  s,  but  the  algorithm  was 
unsuccessful  in  extending  the  curve  to  shorter  peri¬ 
ods.  We  therefore  accumulated  two  data  sets:  a  large 
collection  of  measurements  valid  for  periods  longer 
than  45  s  and  a  slightly  smaller  data  set  extending  to 
the  shortest  period  we  considered,  32  s.  The  quality 
of  each  dispersion  curve  was  assessed  in  various  ways, 
and  the  value  of  the  misfit  function  was  found  to 
provide  th(*  best  indicator  of  overall  quality.  Measure¬ 
ments  with  4^  >0.75  were  found  to  be  less  reliable  and 
were  discarded.  This  led  to  a  reduction  in  the  number 
of  useful  paths  to  approximately  56.000. 

In  further  analysis  of  these  measurements,  the  dis¬ 
persion  curves  were  converted  to  phase  anomalies  5^ 
at  the  di.sciete  periods  150.  100,  75,  60,  50,  45,  40,  37, 
and  35  s.  The  maximum  and  minimum  periods  cor¬ 
respond  roughly  to  the  center  frequencies  of  the  first 
and  hist  misfit  bands  in  equation  (16).  Owing  to  the 
sj)line  parameterization,  the  phase  anomalies  at  ad¬ 
jacent  frequencies  are  correlated.  However,  by  over- 
.sarnpling  the  dispersion  curve  and  calculating  phase 
velocity  maps  at  many  frequencies,  it  will  in  the  future 
be  easier  to  interpolate  smooth  dispersion  curves  for 
arbitrary  paths.  In  the  conversion  to  phase  anomaly 
measurements,  an  additional  quality  check  was  ap¬ 
plied  by  discarding  measurements  with  poor  fit  in  the 
frequency  bands  bracketing  the  desired  frequency. 

The  result  of  applying  these  quality  constraints  is  a 
homogeneous  data  set  of  phase  anomalies.  The  num¬ 
ber  of  observations  varies  depending  on  the  wave  type 
and  period  and  ranges  between  approximately  15,000 
and  38,000  (Table  1).  The  jump  in  the  number  of 


observations  between  the  45  and  50  s  measurements 
is  due  to  the  fact  that  shorter-period  observations  are 
derived  from  the  smaller  data  set  of  dispersion  curves 
described  above.  The  smaller  number  of  good  Love 
wave  observations  at  150  s  reflects  the  difficulty  of 
isolating  the  long-period  fundamental  mode  signal  in 
minor  arc  measurements  at  these  periods. 

In  addition  to  separating  “good”  measurements 
from  “bad”  by  imposing  strict  requirements  on  the 
misfit  function  ^  and  the  fit  in  individual  frequency 
bands,  we  attempted  to  make  realistic  estimates  of 
the  observational  uncertainties  in  our  data  set.  Errors 
in  single-station  phase  velocity  measurements  come 
from  several  sources.  For  example,  if  an  earthquake 
epicenter  is  poorly  known,  an  error  SX  in  the  epi- 
central  distance  will  lead  to  a  phase  error  of  u)5Xlc^, 
which,  for  a  40  km  mislocation  and  a  phase  velocity 
of  4  km  s“^,  leads  to  a  10  s  difference  in  travel  time 
and  to  a  7r/2  phase  error  at  40  s  period.  Errors  in  the 
source  phase  calculation,  due  to  an  incorrect  moment 
tensor  or  earthquake  depth  or  more  subtly  as  a  con¬ 
sequence  of  the  source  excitation  being  calculated  in 
PREM  rather  than  in  a  local  model,  can  lead  to  phase 
errors  as  large  as  tt.  Errors  in  station  timing  are  also 
possible,  as  are  methodological  errors  due  to  inade¬ 
quacies  of  our  dispersion  measurement  algorithm. 

Rather  than  attempting  to  estimate  the  probable 
error  distributions  for  these  different  effects,  we  make 
an  empirical  estimate  of  the  aggregate  uncertainty  in 
our  measurements.  Owing  to  the  scope  of  our  study, 
a  large  number  of  phase  observations  correspond  to 
pairwise  similar  paths.  This  results  primarily  from 
having  many  earthquakes  in  one  region  recorded  at 
an  individual  station,  but  there  are  also  some  clusters 
of  stations,  as  in  California,  which  record  essentially 
the  same  path  for  a  given  earthquake. 

We  first  assume  that  each  phase  anomaly  measure¬ 
ment  for  a  given  path  has  associated  with  it  an 
observational  error  ti  w^hich  is  drawn  from  a  normal 
distribution  with  zero  mean  and  standard  deviation  cr. 
W’e  next  construct  a  distribution  of  measurement  dif¬ 
ferences  by  collecting  all  pairs  of  observations  which 
correspond  to  similar  paths  and  calculate  the  differ¬ 
ence  between  the  two  phase  measurements.  If,  in  each 
pair,  the  two  measurements  are  uncorrelated,  the  re¬ 
sulting  distribution  will  have  a  standard  deviation  of 
2cr.  By  calculating  the  standard  deviation  of  this  dis¬ 
tribution,  w^e  therefore  obtain  an  indirect  measure  of 
2  times  the  uncertainty  in  our  observations.  We  note 
that  the  uncertainty  we  estimate  in  this  way  is  likely 
to  be  smaller  than  the  true  uncertainty,  since  it  does 
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not  include  the  effects  of  systematic  errors,  for  exam¬ 
ple,  in  earthquake  locations. 

For  each  period  that  we  analyzed,  we  constructed 
a  distribution  of  measurement  differences,  selecting 
all  pairs  of  paths  for  which  the  starting  and  ending 
points  w^ere  within  3°  of  each  other.  The  standard  de¬ 
viation  of  this  distribution  was  calculated,  and  half  of 
this  value  was  assigned  as  the  estimated  observational 
uncertainty  in  our  measurements  Figure  6  shows  the 
distribution  of  differences  for  Rayleigh  waves  at  75  s 
period.  On  the  basis  of  our  observation  that  paths 
with  low  misfits  were  of  a  higher  quality,  we  sep¬ 
arated  the  paths  into  three  ranges  of  and  thereby 
obtained  three  groups  of  observations  with  varying 
quality:  The  A  range  (0  <  ^  <  0.25),  the  B  range 
(0.25  <  ^  <  0.5),  and  the  C  range  (0.5  <  <  0.75), 

each  associated  with  an  empirical  uncertainty  (Ta¬ 
ble  2). 

Inversion  for  Global  Phase  Velocity 
Maps 

We  use  the  collection  of  dispersion  measurements 
and  the  associated  uncertainties  to  derive  and  com¬ 
pare  three  models  of  global  surface  wave  propaga¬ 
tion:  a  simple  regionalized  model;  a  low-degree  spher¬ 
ical  harmonic  model  obtained  by  undamped  inversion: 
and  a  smooth  high-degree  spherical  harmonic  model 
obtained  by  damped  inversion.  We  make  use  of  the 
standard  approximation  that  the  surface  waves  have 
traveled  the  minor  arc  along  the  great  circle  connect¬ 
ing  the  source  and  receiver  and  do  not  here  explore 
the  possible  effects  of  more  complicated  propagation 
paths  on  a  heterogeneous  Earth.  Earlier  studies  have 
concluded  that  off-great  circle  propagation  effects  on 
the  propagation  plnuse  usually  are  small,  while  effects 
on  takeoff  and  arrival  azimuths  and  amplitudes  can 
be  significant  [e.g.,  Laske,  1995]. 

A  Pure-Path  Regionalized  Inversion 

In  order  to  quantify  the  quality  of  fit  to  our  data 
that  can  be  obtained  by  a  very  simple  regionalized 
model  of  the  Earth,  we  use  the  dispersion  measure¬ 
ments  to  estimate  phase  velocities  in  the  six  tec¬ 
tonic  regions  identified  by  Jordan  [1981]  in  the  model 
GTRl.  This  geographical  regionalization  provides  a 
gross  differentiation  between  young  (A),  intermedi¬ 
ate  age  (B),  and  old  oceanic  lithosphere  (C),  as  well 
as  between  Precambrian  shields  and  platforms  (S), 
Phanerozoic  platforms  (P),  and  regions  of  the  conti¬ 
nents  having  been  exposed  to  Phanerozoic  deforma¬ 


tion  and  magmatic  activity  (Q),  where  the  letters  cor¬ 
respond  to  the  labels  used  by  Jordan  [1981]. 

We  use  the  pure-path  approach  to  derive  regional 
dispersion  curves.  The  observed  phase  anomaly  is 
modeled  as  having  accumulated  as  a  result  of  phase 
velocity  variations  in  the  six  tectonic  regions, 


j=i 


cP 


(17) 


where  JA*  is  the  path  length  and  is  the  phase  veloc¬ 
ity  variation  in  the  zth  tectonic  region.  We  have  made 
use  of  the  standard  approximation  1/(1  +  Sc/cP)  « 
1  -  Sc/P,  which  is  adequate  as  long  as  the  velocity 
variations  are  small.  In  order  to  find  the  phase  veloc¬ 
ity  variations  5c^  which  optimally  fit  the  observations, 
we  form  the  misfit  function 


where  j  is  the  index  of  the  observation,  N  is  the  to¬ 
tal  number  of  observations,  and  aj  is  the  associated 
observational  uncertainty. 

The  misfit  function  is  minimized  by  straight¬ 
forward  least  squares  inversion  for  the  six  velocity 
perturbations  Jc*  with  respect  to  PREM.  Figure  7 
shows  the  dispersion  curves  for  Love  and  Rayleigh 
waves.  The  results  are  similar  to  those  derived  in  ear¬ 
lier  studies,  and  here  we  focus  on  how  well  the  derived 
regionalized  model  explains  our  observations.  Fig¬ 
ure  8  shows  the  variance  reduction  at  different  peri¬ 
ods.  The  dispersion  curves  derived  for  GTRl  explain 
more  than  80%  of  the  Love  wave  signal  at  the  short¬ 
est  periods  and  approximately  70%  of  the  Rayleigh 
wave  signal.  This  is  primarily  a  consequence  of  the 
strong  signal  generated  by  the  difference  between  con¬ 
tinents  and  oceans.  At  longer  periods,  the  variance 
reduction  diminishes,  and  at  150  s  period,  less  than 
half  of  the  data  variance  for  both  Love  and  Rayleigh 
waves  is  explained  by  GTRl.  This  suggests  that  even 
though  a  strong  correlation  exists,  for  example,  be¬ 
tween  shields  and  old  oceanic  lithosphere  and  fast 
phase  velocities,  this  correlation  is  not  sufficient  to 
explain  the  patterns  of  anomalies.  Figure  9  shows 
the  quantity  yp  IN  ^  which  provides  a  meeisure  of  how 
well  the  model  fits  the  observations,  given  the  esti¬ 
mated  uncertainties  ai  (Table  2).  For  a  good  model, 
X^/A  should  approach  unity.  While  GTRl  explains  a 
significant  fraction  of  the  observed  variance  in  phase 
velocity,  the  GTRl  residuals  are  on  average  \/5  times 
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greater  (for  Rayleigh  waves)  and  \/8  times  greater 
(for  Love  waves)  than  the  a  priori  observational  un¬ 
certainties. 


Models  Expressed  in  Spherical  Harmonics 

We  now  consider  a  smooth,  global  parameteri¬ 
zation  in  terms  of  spherical  harmonics.  The  pre¬ 
dicted  propagation  phase  anomaly  between  a  source 
at  (Os,(t>s)  and  a  receiver  at  (0/2,0i?)  is  written 


6^ 


u!  (5c(0,0) 


=  / 


{Ss^<l>s) 


cO 


ds, 


(19) 


and  the  phase  velocity  perturbation  is  expressed  in 
complex  spherical  harmonics 


dc{0. 0) 


L  m=l 

=  EE 

1=0  m--l 


(20) 


where  the  Yfm  are  the  fully  normalized  surface  spheri¬ 
cal  harmonics  [Edmonds,  1960]  and  L  is  the  maximum 
degree  of  the  expcuision.  The  misfit  function  is  for¬ 
mulated  analogously  to  equation  (18),  and  we  invert 
for  the  coefficients  C/m* 

The  resolution  afforded  by  our  large  set  of  obser¬ 
vations  allows  us  to  invert  for  models  truncated  at 
angular  degree  L  =  16  without  applying  additional 
smoothing  or  regularizations  to  the  solution.  We  use 
the  results  of  this  undamped,  truncated  inversion  in 
two  ways.  First,  the  resulting  models  are  used  to 
calculate  global  rms  velocity  perturbation  values  at 
each  period,  Sc^msIc^K  which  we  then  use  to  quantify 
smoothness  in  our  inversion  for  higher  degree  struc¬ 
ture.  W  e  write 

1/2 

(21) 

which,  with  Aim  =  Ke(C/m)  and  Bim  =  3m(C/m)^ 
becomes 


<5Crms 


1/2 

(22) 


Second,  we  use  the  L  =  16  model  to  identify  out¬ 
liers,  that  is,  observations  which  are  so  poorly  pre¬ 
dicted  by  the  L  =  16  model  that  we  believe  that 
they  are  erroneous.  Fewer  than  2%  of  the  data  were 
identified  as  outliers;  these  data  were  not  used  in  the 
higher-degree  inversions.  However,  for  consistency, 
we  continue  to  include  these  data  in  the  calculations 
of  misfit  and  variance  reduction. 


Not  surprisingly,  the  L  =  16  inversions  provide  a 
much  higher  variance  reduction  and  a  better  fit  to 
the  data  than  does  GTRl  (Figure  8).  Especially  for 
longer  periods,  the  improvement  in  variance  reduction 
is  very  clear,  reflecting  the  lower  correlation  of  mantle 
velocity  anomalies  with  surface  geology.  The  improve¬ 
ment  in  fit  to  the  data,  expressed  as  the  reduction 
in  size  of  the  average  residual,  is  approximately  \/3 
(Figure  9).  Figure  10  shows  the  resulting  global  rms 
velocity  variations  for  the  range  of  periods  considered. 
The  effect  of  the  crust  is  seen  in  the  rapid  increase  of 
Love  wave  rms  velocity  variations  at  shorter  periods. 
At  75  s.  Love  and  Rayleigh  waves  exhibit  approxi¬ 
mately  equal  lateral  variability. 

Regularization  of  the  inversion  is  required  when  we 
determine  phase  velocity  maps  expanded  to  higher 
angular  degrees.  Several  choices  are  available  for  sta¬ 
bilizing  the  inversion,  and  we  have  chosen  to  mini¬ 
mize  a  measure  of  the  model  roughness,  defined  as 
the  squared  rms  gradient  of  the  model.  The  scaled 
rms  gradient  is 


(23) 


and  for  our  models  expanded  in  spherical  harmonics, 
we  calculate  TZ  as 


n  = 


(fJCrms/c®) 


l=l 


E  2-4L  + 


(24) 


The  effect  of  smoothing  is  similar  to  applying  a  low- 
pass  filter  to  the  model,  and  an  alternative  defini¬ 
tion  of  roughness,  such  as  the  squared  Laplacian  [e.g., 
Laske  and  Masters,  1996],  corresponds  to  choosing  a 
different  spectral  falloff. 

We  determine  the  model  coefficients  which  solve 
the  minimization  problem 


(25) 


where  7  is  a  damping  parameter  which  needs  to  be 
chosen.  Through  scaling  of  the  gradient  by  {Scrms/c^)~^  5 
we  compensate  for  the  fact  that  shorter-period  maps, 
because  of  their  greater  amplitude  variations,  also 
have  larger  gradients.  This  scaling  allows  us  to  use 
the  same  damping  parameter  7  at  all  periods.  The 
choice  of  7  is  usually  based  on  an  examination  of  the 
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trade-off  curve  between  misfit  and  either  the  size  or 
the  roughness  of  the  resulting  model.  Our  choice  is 
based  on  these  considerations,  aided  by  our  subjec¬ 
tive  judgment  of  what  constitutes  a  reasonable  level 
of  roughness.  Figure  11  shows  a  representative  trade¬ 
off  curve  and  indicates  the  solution  which  corresponds 
to  our  final  choice  for  7. 

Plates  la  and  lb  show  the  resulting  Love  and 
Rayleigh  wave  phase  velocity  perturbation  maps  at 
35,  50,  75,  and  150  s  period.  Many  of  the  features  ap¬ 
parent  in  these  maps  have  been  seen  in  earlier  stud¬ 
ies  as  well  and  have  been  related  to  global  tectonic 
features.  For  example,  there  are  strong  signals  asso¬ 
ciated  with,  at  shorter  periods,  the  contrast  between 
continents  and  oceans  and,  at  all  periods,  slow  ve¬ 
locities  near  Iceland  and  the  Red  Sea  rift  zone.  The 
expected  pattern  of  increasing  velocities  in  the  ocean 
basins,  related  to  the  cooling  of  the  oceanic  litho¬ 
sphere,  is  also  clearly  apparent.  Our  maps  offer  a 
sharper  definition  of  many  features,  allowing  for  a 
firmer  association  between  observed  anomalies  and 
specific  tectonic  elements.  For  example,  for  Love 
waves  at  150  s  and  Rayleigh  waves  at  75  s,  w^e  see 
a  strong  correlation  between  fast  velocities  and  the 
oldest  cratons  and  the  distinct  separation  of  these 
anomalies  along  Proterozoic  orogens;  note  the  separa¬ 
tion  of  the  West  Africa,  Congo,  and  Kalahari  cratons. 
Another  interesting  feature  of  the  maps  is  the  slow 
anomalies  in  the  Pacific  Ocean,  away  from  the  ridges, 
seen  most  prominently  at  the  longest  periods.  A  de¬ 
tailed  interpretation  of  the  many  intriguing  features 
identifiable  in  these  maps  of  phase  velocity,  particu¬ 
larly  in  terms  of  a  3-D  velocity  model  for  the  upper 
mantle,  will  be  presented  in  a  separate  paper. 

Table  3  gives  the  integral  measures  of  rms  level  of 
the  models  and  the  model  gradients  for  the  L  =  40,  . 
the  L  =  16,  and  other  models,  allowing  for  a  compar¬ 
ison  of  their  amplitude  and  roughness.  The  correla¬ 
tions  between  the  L  =  40  and  L  =  16  maps  are  very 
high  and  vary  from  0.95  to  0.99  depending  on  the  pe¬ 
riod  and  wave  type  (Table  4).  These  high  values  are. 
to  a  large  extent,  due  to  the  fact  that  the  heterogene¬ 
ity  in  the  maps  is  dominated  by  low-degree  structures. 
However,  correlations  for  individual  degrees  are  also 
high,  ranging  from  0.67-1.00.  For  0<  /  <  8,  the  cor¬ 
relations  are  all  greater  than  0.95. 

The  model  coefficients  are  not  tabulated  here,  but 
they  can  be  obtained  directly  from  the  authors  or  in 
electronic  form  via  anonymous  ftp  at  saf  .  harvard .  edu, 
by  retrieving  the  files  in  pub/ETL-JGR-96.  The  tabu¬ 
lated  values  are  -4/^  and  Bim  for  m  >  0. 


A  Resolution  Test 

It  is  evident  from  Figures  8  and  9  that  the  im¬ 
provement  in  variance  reduction  and  data  fit  that 
correspond  to  going  from  the  L  =  16  (289  parame¬ 
ters)  models  to  the  L  =  40  (1681  parameters)  mod¬ 
els  is  modest,  and  it  is  therefore  necessary  to  justify 
the  added  model  complexity.  We  performed  several 
standard  resolution  experiments,  which  showed  that 
the  geometry  and  amplitude  of  low-degree  (/  <  12) 
structure  is  essentially  completely  recovered  in  the 
inversions.  Higher-degree  structure  is  recovered  bet¬ 
ter  in  the  northern  than  in  the  southern  hemisphere, 
reflecting  the  ray  coverage.  The  geometry  of  higher- 
degree  structure  is  generally  preserved,  but  the  damp¬ 
ing  used  in  the  inversion,  as  well  as  the  intrinsic  path¬ 
averaging  effect  of  surface  wave  dispersion  observa¬ 
tions,  lead  to  a  smearing  and  reduction  in  amplitude 
of  shorter- wavelength  features.  These  general  conclu¬ 
sions  are  well  illustrated  by  the  following  test  in  which 
w'e  attempted  to  resolve  a  realistic  phase  velocity  map 
based  on  an  a  priori  model  of  the  crust. 

The  thickness  and  velocity  structure  of  the  crust 
contribute  significantly  to  the  dispersion  of  surface 
waves,  in  particular,  at  shorter  periods.  Recently, 
Mooney  et  ai  [1995]  compiled  a  detailed  preliminary 
global  crustal  structure  model,  CRUST-5.0,  with  5° 
by  5°  resolution.  We  calculated  the  dispersion  pre¬ 
dicted  by  this  model  by  substituting  the  CRUST-5.0 
values  {vp,  vs,  and  p)  into  PREM,  adjusting  also  the 
depth  of  the  Moho  and  accounting  for  topography  and 
bathymetry.  We  then  evaluated  the  surface  wave  dis¬ 
persion  in  each  modified  local  Earth  structure.  The 
5°  by  5°  global  dispersion  map  that  resulted  wets  then 
expanded  in  spherical  harmonics  up  to  degree  /  =  40. 
Figure  12  shows  the  phase  velocity  map  thus  calcu¬ 
lated  for  Love  weaves  at  35  s  period. 

For  the  collection  of  paths  that  were  used  in  our  in¬ 
versions,  we  next  calculated  the  corresponding  phase 
anomalies  for  propagation  in  the  synthetic  model.  To 
each  observation,  we  added  Gaussian  noise,  corre¬ 
sponding  to  the  uncertainty  ai  associated  with  each 
path,  multiplied  by  a  factor  of  2  in  order  to  mimic  the 
quality  of  fit  obtained  using  the  real  data.  We  then 
inverted  the  simulated  data  using  the  same  damp¬ 
ing  7  as  w'as  used  in  the  inversion  of  the  real  obser¬ 
vations.  We  achieved  a  96%  variance  reduction  and 

/N=l. Si,  similar  to  what  we  obtained  for  the  real 
data  (Figures  8  and  9). 

Figure  12  shows  the  inversion  results  for  Love 
waves  at  35  s  period.  There  is  very  good  agreement 
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between  the  input  and  retrieved  maps,  in  particular, 
in  the  northern  hemisphere.  An  area  which  shows 
clear  indications  of  smearing  is  the  Andes  mountains, 
w^here  the  narrow’  band  of  north-south  trending  slow' 
velocities  is  poorly  recovered.  In  contrast,  slow  veloc¬ 
ities  corresponding  to  areas  of  thick  crust  and  deep 
basins  in  Eurasia  are  very  well  recovered.  The  essen¬ 
tially  uniform  velocity  of  the  ocean  basins  is  reason¬ 
ably  well  recovered,  with  spurious  variations  not  ex¬ 
ceeding  ±1.5%  (corresponding  to  approximately  10% 
of  the  full  range  of  variations  seen  in  the  map).  Fig¬ 
ure  13  show's  a  comparison  of  the  spectral  pow’ers  of 
the  input  and  retrieved  models,  as  well  as  the  corre¬ 
lation  coefficient  for  each  spherical  harmonic  degree. 
The  power  in  the  input  model  is  well  recovered  up  to 
/  =  12.  and  for  higher  angular  degrees,  the  retrieved 
rno(l(‘l  underestimates  the  pow’er  by  a  factor  as  high  as 
3.  Tin*  corr(‘lation  betw'een  input  and  retrieved  mod¬ 
els  is  greater  than  0.5  for  I  <36.  Since  we  have  fewer 
obs(*r vat  ions  for  Love  w’aves  at  35  s  than  for  any  other 
inversion,  we  believe  that  this  test  provides  a  good  in¬ 
dication  of  how  well  w  e  recover  the  true  phase  velocity 
variations  in  the  Earth. 

Discussion 

Plate  1  and  the  corresponding  spherical  harmonic 
coefficients  repr(*sent  the  primary  results  of  our  study. 
Interi)retation  of  these  results  in  terms  of  3-D  Earth 
structure  will  b(‘  presented  elsew’here.  Here  w*e  limit 
our  discu.ssion  to  (1)  the  discrepancy  between  our 
global,  average  pha.se  velocities  and  those  of  PREM. 
(2)  a  comparison  of  our  maps  with  the  effects  of 
crustal  structure,  and  (3)  comparisons  w'ith  other  re¬ 
cent  studies  of  int(*r mediate- period  global  plnise  ve¬ 
locity  variations. 

Spherically  Symmetric  Term 

The  results  of  our  study  indicate  that  the  spluTic  al 
average  phase  velocities  predicted  by  PREM  require* 
large  corrections.  Figure  14  show’s  the  spherical  av¬ 
erage  perturbation  for  Love  and  Rayleigh  waves  at 
different  periods,  calculated  from  the  Aqo  term  in  the* 
spherical  harmonic  expansions.  For  Rayleigh  wavers 
around  50  s,  the  correction  required  is  about  0.8% 
w’hich  is  nearly  half  of  the  observed  rms  variation 
(Figure  14).  We  note  that  the  sign  of  the  correction 
is  the  same  as  that  observed  betw'een  data  used  to  dep¬ 
rive  PREM  and  PREM  predictions  [Dziewon.ski  and 
Anderson,  1981].  However,  using  these  new  data,  the 
magnitude  of  the  difference  is  greater.  While  there 


are  many  possible  explanations  for  the  observed  dis¬ 
crepancy,  we  suggest  that  since  the  Love  wave  data 
are  fit  quite  well  by  PREM  and  the  Rayleigh  waves 
are  significantly  faster  than  predicted,  the  transverse 
anisotropy  at  the  top  of  the  mantle  may  not  be  as 
strong  as  in  PREM,  where  it  reaches  5%  just  below’ 
the  Moho.  Other  explanations  related,  for  example, 
to  the  global  average  crustal  structure,  are  of  course 
also  plausible.  In  this  article,  we  have  not  included 
the  spherically  symmetric  term  in  the  calculations  of 
rms  values,  variance  reductions,  or  correlations  be¬ 
tween  models. 

Effect  of  the  Crust 

In  tomographic  studies  of  the  Earth’s  mantle,  sur¬ 
face  wave  dispersion  mecisurements  and  phase  velocity 
maps  are  commonly  used  to  constrain  upper  mantle 
velocity  variations.  For  this  use,  it  is  necessary  to 
account  for  crustal  effects  on  the  observations  before 
inferring  from  them  the  nature  of  the  deeper  struc¬ 
ture,  This  is  commonly  done  by  subtracting  out  the 
predicted  effects  of  an  a  priori  crustal  structure,  often 
one  W’hich  corresponds  to  the  average  oceanic  or  con¬ 
tinental  crust.  The  crustal  signal  is  strong  for  shorter- 
period  w'aves  and  decreases  at  longer  periods.  Since 
the  rms  variations  in  observed  phase  velocities  also  de- 
crea.se  with  increasing  period,  it  is  not  clear,  however, 
that  the  relative  size  of  the  crustal  signal  decreases 
with  increasing  period.  We  use  the  model  predictions 
bascM:!  on  Mooney  et  a/.’s  [1995]  crustal  model,  de¬ 
scribed  above,  to  investigate  this  relationship. 

The  correlation  betw’een  observed  phase  velocities 
and  those  due  to  the  crust  (Figure  15a)  follows  the 
expected  trend  of  a  greater  effect  at  short  periods 
*md  for  Love  waves.  The  0.83  correlation  for  Love 
waves  at  35  s  show's  that  while  Love  w’aves,  even  at 
this  period,  have  some  sensitivity  to  mantle  structure, 
the  observed  variations  are  dominated  by  the  crust. 
The  good  agreement  in  the  patterns  of  the  “observed” 
and  “predicted”  maps  is  also  seen  in  the  amplitudes 
of  their  rms  variations;  the  rms  variation  (Figure  15b) 
in  the  map  predicted  from  CRUST-5.0  is  4.48%  and 
from  the  inverted  phase  velocity  maps  is  3,91%.  In 
addition,  the  rms  gradient  of  our  model  TZ  =  78.5  is 
similar  to  the  7^  =  91.0  we  find  for  the  crustal  model. 

The  correlation  with  the  crustal  model  decreases 
with  increasing  period,  and  for  periods  longer  than 
50  s  for  Rayleigh  waves  and  100  s  for  Love  waves,  the 
correlation  is  negative  (Figure  15a).  When  we  cal¬ 
culate  the  correlation  for  continental  regions  only,  it 
is  seen  to  deteriorate  faster  with  increasing  period. 
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so  that  for  Rayleigh  waves,  the  correlation  becomes 
negative  at  around  45  s,  presumably  reflecting  a  lack 
of  correlation  between  the  thickness  of  the  continen¬ 
tal  crust  and  the  seismic  velocities  of  the  underlying 
mantle. 

Figures  15b  and  15c  compare  the  rms  levels  of  our 
phase  velocity  maps  with  those  of  the  crustal  model 
riding  on  a  spherically  symmetric  mantle.  Also  shown 
is  the  '  ms  variation  of  the  residual  model  which  is  ob¬ 
tained  if  the  crustal  phase  velocity  maps  (as  derived 
from  Mooney  et  a/.’s  [1995]  model)  are  subtracted 
from  our  observed  maps.  The  level  of  heterogeneity 
seen  for  Love  waves  is  similar  to  that  which  we  would 
expect  from  the  crust  alone.  Of  course,  because  of  the 
lack  of  correlation  at  longer  periods,  there  is  a  mantle 
signal  in  the  observed  maps,  but  it  is  not  dominant. 
For  Rayleigh  waves,  the  observed  phase  velocity  varia¬ 
tions  are  greater  than  those  due  only  to  the  crust,  but 
it  is  interesting  to  note  that  the  level  of  variations  is 
most  similar  at  the  limits  of  our  period  interval.  This 
suggests  that  the  details  of  the  crustal  correction  that 
are  applied  in  tomographic  inversions  for  3-D  struc¬ 
ture  are  important  not  only  at  short  periods  but  also 
at  100  s  and  greater  periods. 

Comparison  With  Other  Studies 

Several  groups  are  pursuing  research  similar  to  that 
presented  here,  and  it  useful  to  compare  our  results 
with  other  published  studies.  At  periods  longer  than 
80  s,  there  exists  a  long  history  of  global  phase 
locity  models,  but  since  many  of  the  earlier  stud¬ 
ies  have  been  superseded  by  new  and  more  complete 
studies,  we  limit  the  comparisons  to  four  recent  stud¬ 
ies:  Trampert  and  Woodhouse  [1995,  1996],  Zhang  and 
Lay  [1996],  and  Laske  and  Masters  [1996].  In  particu¬ 
lar,  the  two  studies  by  Trampert  and  Woodhou.se  are 
similar  in  scope  to  the  work  presented  here.  Tram- 
pert  and  Woodhouse  [1995]  applied  an  automatic  al¬ 
gorithm  to  analyze  the  dispersion  of  a  large  number  of 
surface  waves  and  derived  global  phase  velocity  maps 
expanded  up  to  spherical  harmonic  degree  /  =  40. 
A  difference  in  the  data  used  between  our  studies  is 
that  Trampert  and  Woodhouse  [1995,  1996]  used  both 
minor  and  major  arc  arrivals,  while  we  use  only  the 
minor  arc.  We  derive  results  to  35  s  period,  while 
Trampert  and  Woodhouse's  [1995,  1996]  shortest  p(^ 
riod  is  40  s.  In  the  inversion  for  phase  velocity  maps, 
Trampert  and  Woodhouse  regularized  the  solution  by 
minimizing  the  squared  Laplacian,  in  contrast  to  our 
choice  of  minimizing  the  squared  gradient.  In  their 
two  papers,  in  1995  and  1996,  results  are  presented 


for  40,  60,  80,  100,  and  150  s  periods  for  Love  and 
Rayleigh  waves.  The  primary  differences  between  the 
two  Trampert  and  Woodhouse  studies  are  the  qual¬ 
ity  and  quantity  of  data  included  in  the  inversions 
and  the  level  of  damping.  The  1996  results  are  less 
damped  and  consequently  rougher.  Since  our  results 
are  in  better  agreement  with  the  earlier  results,  we 
make  comparisons  with  both  studies.  Zhang  and  Lay 
[1996]  derived  phase  velocity  maps  between  85  and 
250  s  using  both  minor  and  major  arc  data.  We  com¬ 
pare  our  results  with  theirs  at  100  s  and  150  s.  Laske 
and  Masters  [1996]  similarly  derived  maps  between 
~75  and  250  s.  We  compare  our  maps  with  theirs  at 
75,  100,  and  150  s. 

Table  3  gives  the  rms  and  rms  gradient  values  for 
the  various  maps,  and  Table  4  summarizes  their  cor¬ 
relations  with  our  preferred  L  =  40  maps.  For  Love 
waves  at  100  and  150  s,  the  three  maps  by  Zhang  and 
Lay  [1996],  Laske  and  Masters  [1996],  and  Trampert 
and  Woodhouse  [1995]  correlate  equally  well  with  our 
maps.  No  two  maps  from  these  three  studies  corre¬ 
late  better  than  they  do  individually  with  our  map. 
For  Rayleigh  waves,  our  results  have  the  best  cor¬ 
relation  with  Laske  and  Masters’s  results.  As  can 
be  seen,  we  have  consistently  better  agreement  with 
Trampert  and  Woodhouse’s  1995  study  than  with  the 
1996  study. 

The  correlation  of  Trampert  and  Woodhouse'^s  [1996] 
Love  wave  map  at  40  s  period  with  ours  is  0.80,  sug¬ 
gesting  a  relatively  good  agreement  between  the  two 
models.  However,  it  is  worth  noting  that  we  ob¬ 
tained  an  almost  equally  good  correlation  (0.79)  be¬ 
tween  our  map  and  that  predicted  from  the  a  priori 
crustal  model  (see  Figure  15).  Thus  the  correlation 
between  the  two  results  is  primarily  due  to  the  large- 
scale  spatial  pattern  of  oceanic  and  continental  crust, 
and  there  exist  significant  disagreements  between  the 
two  maps  concerning  smaller  structures.  Figure  16 
shows  the  correlation  between  the  two  models  for  dif- 
f<‘r(*nt  angular  degrees  in  the  spherical  harmonic  ex¬ 
pansion.  The  correlations  between  the  two  models 
start  to  diverge  for  degrees  as  low  as  /  =  7,  and  while 
they  remain  positive  for  most  angular  degrees,  they 
approach  zero  for  the  highest  degrees. 

We  believe  that,  at  shorter  periods,  our  results  of¬ 
fer  an  improvement  over  the  models  of  Trampert  and 
Woodhouse  [1995,  1996],  based  on  two  observations. 
First,  the  reductions  in  variance  reported  by  Tram¬ 
pert  and  Woodhouse  are  low  in  comparison  with  what 
we  obtain  in  our  study.  For  example,  we  obtain  a  vari¬ 
ance  reduction  of  96%  of  our  Love  wave  observations 
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at  40  s  period,  while  Trampert  and  Woodhouse  in 
their  1995  and  1996  studies  report  86%  and  77%,  re¬ 
spectively.  In  fact,  the  fits  that  their  models  provide 
to  our  data  are  better  (90%  and  87%,  respectively) 
than  their  reported  values.  This,  to  us,  provides  a 
strong  suggestion  that  their  raw  measurements  are 
significantly  noisier  than  ours  and  internally  less  con¬ 
sistent,  perhaps  suggesting  a  problem  with  the  mea¬ 
surement  technique. 

The  second  observation  we  make  is  that  while  in 
many  regions  there  is  a  striking  similarity  between 
our  two  models,  there  are  some  areas  where  there  is 
strong  disagreement,  and  we  believe  our  result  is  more 
realistic.  Figure  17  shows  a  comparison  between  two 
Love  wave  maps  at  40  s  for  India  and  central  Asia,  a 
region  which  contains  the  largest  elevated  plateau  and 
tlie  thickest  crust  on  Earth.  These  geological  features 
are  reflected  in  our  map  by  the  extremely  slow'  phase 
velocity  across  Tibet  and  the  continuous  band  of  slow 
velocities  along  the  Alpine-Himalayan  belt  of  thick¬ 
ened  and  tectonically  active  crust.  In  Trampert  and 
Woodhouse's  [1996]  map  and  similarly  in  their  1995 
study,  the  Tibetan  plateau  is  not  as  prominent  a  fea¬ 
ture,  and  it  is  cut  and  truncated  by  a  fast  anomaly  at 
approximately  85°E.  On  the  basis  of  what  is  knowm 
about  the  cru.st  and  tectonics  for  this  region  (see  Fig¬ 
ure  12).  as  well  as  the  results  from  tomographic  stud¬ 
ies  [e.g..  Bourjot  and  Romanowicz,  1992],  it  appears 
highly  unlikely  that  Trampert  and  W’oodhouse’s  fast 
anomaly  is  a  real  feature  of  the  Earth.  The  fact  that 
southern  Eurasia  is  quite  well  sampled  by  crossing 
wave  paths  in  both  studies  also  suggests  that  the  dif¬ 
ferences  are  not  due  to  limited  resolution  but  are  in¬ 
stead  due  to  errors  in  the  measuren?ents. 

Conclusions 

We  have  developed  a  new  algorithm  for  making 
automated  phase  measurements  of  dispersed  surface 
waves.  By  insisting  that  the  surface  w'ave  signal 
should  have  a  smoothly  varying  phase  and  amplitude 
and  by  anchoring  the  dispersion  curve  at  a  period 
w’here  there  is  no  ambiguity  in  the  total  phase,  the 
new  method  yields  robust  measurements  of  phase  for 
Love  and  Rayleigh  waves  with  periods  as  short  as 
35  s.  It  is  noteworthy  that  some  of  these  phases  ex¬ 
hibit  phase  delays  or  advances  of  as  many  as  10  cycles 
(±350  s).  The  potential  exists  for  applying  this  algo¬ 
rithm  to  higher  frequencies  as  w'ell,  w'hich  might  be 
of  particular  use  for  analysis  of  recordings  at  regional 
distances.  ^ 


The  systematic  application  of  our  algorithm  to 
teleseismic  data  from  the  Global  Seismographic  Net¬ 
work  has  yielded  a  very  large  database  of  dispersion 
and  amplitude  measurements.  By  comparing  results 
corresponding  to  similar  paths,  we  have  been  able  to 
estimate  the  uncertainties  in  our  observations.  These 
uncertainties  are  useful,  since  without  them  it  is  easy 
to  erroneously  conclude  that  a  given  model  provides 
a  good  fit  to  the  data  based  on  a  high  variance  reduc¬ 
tion.  For  example,  one  might  be  led  to  believe  that 
the  simple  regionalization  GTRl  provides  a  good  fit 
to  shorter-period  surface  waves  since  it  reduces  the 
variance  by  65-85%  (Figure  7).  In  fact,  while  GTRl 
explains  a  large  fraction  of  the  signal  in  the  data  at 
short  periods,  primarily  due  to  the  large  contrast  be¬ 
tween  continental  and  oceanic  structures,  the  fit  it 
provides  to  the  observations  is  mediocre  at  all  peri¬ 
ods.  Knowledge  of  the  approximate  true  uncertain¬ 
ties  in  our  data  also  allowed  us  to  assess  the  resolu¬ 
tion  and  recovery  of  true  Earth  structures  through  a 
realistic  synthetic  test  (Figures  12  and  13).  Our  con¬ 
clusions  differ  substantially  from  those  of  Ricard  et  al 
[1996],  who  performed  a  similar  test  using  the  a  priori 
model  3SMAC  [Nataf  and  Ricard^  1996]  and  a  differ¬ 
ent  path  coverage.  Ricard  et  al  [1996]  obtained  a  very 
low'  recovery  of  pow'er  in  degrees  greater  than  10  and 
only  a  53%  recovery  of  the  signal  for  degree  3.  They 
concluded,  in  agreement  with  some  previous  studies 
[Mochizuki,  1993;  Snieder,  1993;  Nolet  et  al,  1994], 
that  the  spectrum  of  lateral  heterogeneity  cannot  be 
deduced  for  degrees  higher  than  10  using  global  sur¬ 
face  wave  tomography.  In  contrast,  we  obtain  good 
recovery  of  structures  (Figure  13)  to  at  least  degree  20 
and  essentially  complete  recovery  of  the  structure  up 
to  degree  12.  While  we  do  underestimate  the  power 
of  higher  order  structures,  this  is  not  as  severe  a  prob¬ 
lem  as  implied  by  the  above  authors.  We  agree,  how¬ 
ever,  that  the  best  way  to  achieve  better  recovery  of 
short er-w'avelength  structures  is  to  include  observa¬ 
tions  corresponding  to  shorter  paths,  provided  the 
difficulties  of  making  these  measurements  are  over¬ 
come.  In  addition,  data  which  are  more  sensitive  to 
gradients  in  the  models,  such  as  polarization  angles 
and  amplitudes,  could  help  improve  the  resolution  of 
smaller-scale  structures  [Laske  and  Masters,  1996]. 

A  central  contribution  of  this  paper  is  the  phase  ve¬ 
locity  maps.  These  were  derived  using  several  simpli¬ 
fying  assumptions.  First,  we  assumed  that  the  phase 
measurements  could  be  interpreted  using  the  unper¬ 
turbed  great  circle  path.  Clearly,  this  assumption  is 
easily  challenged,  given  the  level  of  heterogeneity  re- 
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trieved,  for  example,  for  40  s  Love  waves  in  Tibet 
(Figure  17).  On  the  basis  of  ray-tracing  experiments, 
we  do  not  believe  that  taking  more  accurate  account 
of  the  refraction  caused  by  such  heterogeneity  will  af¬ 
fect  the  larger-scale  features  in  our  models.  It  is  clear, 
however,  that  for  certain  paths,  such  as  those  crossing 
Tibet  or  grazing  continent-ocean  margins,  the  sys¬ 
tematic  bending  of  rays  into  slow  areas  may  affect 
the  mapping  of  boundaries  and  smaller-scale  features. 
Further  research  is  needed  to  investigate  the  poten¬ 
tial  bias  introduced  by  this  simplifying  assumption  in 
specific  areas. 

A  second  important  assumption  is  that  the  ob¬ 
served  dispersion  can  be  explained  by  isotropic  varia¬ 
tions  in  phase  velocity.  Several  previous  studies  have 
concluded  that  azimuthal  anisotropy  contributes  sig¬ 
nificantly  to  the  observed  variability  in  surface  wave 
dispersion,  in  particular,  for  Rayleigh  waves  travers¬ 
ing  oceanic  lithosphere  [e.g.,  Forsyth,  1975;  Tanimoto 
and  Anderson,  1985;  Nishimura  and  Forsyth,  1988: 
Montagner  and  Tanimoto,  1990].  With  the  database 
collected  in  this  study,  we  are  currently  investigat¬ 
ing  the  effect  of  aizimuthal  anisotropy.  Like  previous 
investigators,  we  find  that  the  main  difficulty  lies  in 
distinguishing  between  isotropic  and  anisotropic  het¬ 
erogeneity  in  areas  with  poor  path  coverage.  Results 
from  this  investigation  w’ill  be  presented  elsewhere. 
Here  it  suffices  to  say  that  the  good  fit  to  the  data  that 
can  be  achieved  without  accounting  for  anisotropy 
suggests  that  anisotropy,  while  probably  important 
in  certain  regions,  remains  a  second-order  effect  when 
considering  a  global  data  set. 

Notwithstanding  the  need  to  investigate  further 
the  limitations  of  the  simplifying  assumptions  em¬ 
ployed  in  our  inversions,  the  phase  velocity  maps  de¬ 
rived  here  provide  very  good  fits  to  the  global  data 
set  we  have  collected.  The  maps  can  therefore  also  be 
used  to  predict  the  phase  of  teleseLsmic  surface  waves, 
something  which  will  be  useful  for  earthquake  stud¬ 
ies.  The  most  important  use  for  these  maps  and  the 
dispersion  data  is,  however,  in  the  determination  of 
the  3-D  structure  of  the  upper  mantle.  Preliminary 
analysis  of  our  data  set  in  the  development  of  a  model 
expanded  up  to  degree  20  in  the  upper  mantle  [Ek^ 
Strom  and  Dziewonski,  1995,  1996]  indicates  that  our 
measurements  of  Love  and  Rayleigh  wave  dispersion 
are  consistent  with  previous  data  sets  but  that  they 
require  higher  amplitudes  of  lateral  heterogeneity  in 
the  shallow’est  mantle  than  seen  in  previous  Harvard 
mantle  models  [e.g.,  Su  et  al,  1994]. 

The  experiments  described  above  have  shown  that 


the  crustal  signal  in  the  phase  velocity  maps  is  sig¬ 
nificant  and  that  it  must  be  adequately  removed  be¬ 
fore  inferring  the  structure  of  the  uppermost  mantle. 
While  Mooney  et  a/.’s  [1995]  model  clearly  gives  a  bet¬ 
ter  picture  of  the  global  crust  than  any  we  have  had 
previously,  the  differences  between  Figures  12  and  17 
appear  to  be  difficult  to  explain  without  adjustments 
to  the  crustal  structure.  A  future  inversion  of  our 
data  which  simultaneously  adjusts  both  crustal  and 
mantle  3-D  structure  therefore  seems  most  desirable. 
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Figure  1.  The  six  B-spline  basis  functions  fi  used  to  parameterize  the  apparent  phase  velocity  perturbation  Sc 
and  the  spectral  amplitude  A. 

Figure  2.  Eleven  partially  overlapping  frequency  bands,  evenly  spaced  in  frequency  between  250  and  32  s  period. 

Figure  3.  Comparison  of  observed  and  model  seismograms  at  progressive  stages  of  the  inversion,  (a)  Initial  fit, 
before  adjustment  of  phase  velocity  and  amplitude,  250-32  s.  (b)  Initial  fit,  250-75  s.  While  the  misfit  is  large,  it 
is  apparent  that  the  waveform  can  be  matched  after  minor  adjustments  of  the  phase,  (c)  Fit  after  adjustment  of 
(Sc(cj)  and  -4(a;)  in  the  period  band  250-75  s.  (d)  Fit  after  adjustment  in  the  period  band  250-40  s.  (e)  Final  result 
after  fitting  the  data  in  the  period  band  250-32  s. 

Figure  4.  Dispersion  and  spectral  amplitude  results  for  the  analysis  summarized  in  Figure  3.  (a)  Initial  guess 
(from  SH8U4L8)  and  final  result  for  the  phase  velocity  perturbation  expressed  as  a  percentage  of  the  preliminary 
reference  Earth  model  (PREM)  phase  velocity  at  each  frequency,  (b)  Initial  and  final  amplitude  spectra,  normalized 
to  the  largest  value  in  the  range  of  frequencies. 

Figure  5.  Map  showing  the  locations  of  1744  earthquakes  (hexagons)  and  158  stations  (squares),  each  of  which 
forms  the  start  or  end  point  for  at  least  one  successfully  mecisured  Love  or  Rayleigh  wave  path. 

Figure  6.  Distribution  of  phase  measurement  differences  for  Rayleigh  waves  at  75  s  and  pairwise  similar  paths. 
The  thro(‘  different  curves  correspond  to  measurements  of  A,  B,  2u:id  C  quality.  The  curves  have  been  normalized 
to  have  the  same  area. 

Figure  7.  Dispersion,  expressed  as  perturbations  with  respect  to  PREM,  determined  for  the  six  tectonic  regions 
ofGTRl  [Jordan,  1981].  (a)  Love  waves,  (b)  Rayleigh  waves. 

Figure  8.  \*ariance  reduction  of  the  complete  data  set  and  three  global  models:  the  GTRl  regionalization  and 
dispersion  curves,  the  L  =  16  undamped  spherical  harmonic  model,  and  the  preferred  L  =  iO  smooth  model. 

(a)  Love  wav(‘s.  (b)  Rayleigh  waves. 

Figure  9.  Tlu*  goodness  of  fit  parameter  calculated  for  different  inversions  and  models:  the  GTRl  regional¬ 

ization  and  dispersion  curves,  the  L  =  16  undamped  spherical  harmonic  model,  and  the  preferred  L  =  40  smooth 
model,  (a)  Love*  waves,  (b)  Rayleigh  waves. 

Figure  10.  The  rms  \'alue  of  the  velocity  perturbations  at  different  periods  for  the  undamped  L  =  16  inversions. 
The  /  =  0  term  is  not  included  in  the  rms  value. 

Figure  11.  Trach^off  curve  for  the  L  =  40  inversion  of  Love  waves  at  60  s  period,  showing  the  effect  of  choosing 
different  values  of  in  cHjuation  (25).  The  shape  of  the  curve  is  representative  of  a\V L  =  40  inversions. 

Figure  12.  (top)  The  Love  wave  phase  velocity  at  35  s  predicted  from  Mooney  et  al's  [1995]  model  of  the  crust, 
(bottom)  The  recovery  of  this  model  after  synthesizing  obser\'ations  corresponding  to  the  actual  path  coverage, 
adding  realistic  noise,  and  performing  the  danqxHl  inversion. 

Figure  13.  Graph  showing  the  recovery  of  individual  spherical  harmonic  degrees  in  the  resolution  test,  (a)  Power 
of  individual  degrees  in  the  model.  Solid  triangles  corrt^spond  to  the  input  model;  open  ones  correspond  to  the 
retrieved  model,  (b)  Correlation  between  the  input  and  retrieved  models  as  a  function  of  angular  degree. 

Figure  14.  The  spherical  average  phase  velocity  perturbation  with  respect  to  PREM  for  Love  and  Rayleigh  waves. 

Figure  15.  (a)  Correlation  between  our  preferr(‘d  L  =  40  maps  and  the  phase  velocity  maps  calculated  from  the 
crustal  model  of  Mooney  et  al  [1995].  Both  the  global  value  and  the  value  for  continental  areas  alone  are  shown. 

(b)  The  observed  rms  value  of  Love  wave  pha.se  velocity  variations  compared  with  the  predictions  from  Mooney  et 
alJs  [1995]  crustal  model.  The  rms  value  of  the  difference  between  observed  and  predicted  phase  velocity  variations 
is  also  shown,  (c)  Same  as  (b),  except  for  Rayleigh  waves. 
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Figure  16.  Comparison  of  Trampert  and  Woodhouse  s  [1996]  Love  wave  map  at  40  s  and  ours,  (a)  Correlation  of 
individual  spherical  harmonic  degrees,  (b)  Power  of  the  two  models  at  degrees  1-40. 

Figime  17.  Regional  map  showing  the  Love  wave  phase  velocity  at  40  s  period,  (top)  The  map  of  Trampert  and 
Woodhouse  [1996]  and  (bottom)  our  results.  In  general,  there  is  less  correlation  in  the  top  map  with  known  crustal 
and  tectonic  structures,  such  as  the  thick  crust  of  Tibetan  platform  and  the  continent-ocean  contrast  around  India. 
In  other  areas,  for  example,  east  of  100°E,  the  two  models  are  extremely  similar.  The  overall  correlation  between 
the  two  regional  maps  is  0.68. 


Plate  la.  Global  Love  wave  phase  velocity  maps  at  35,  50,  75,  and  150  s.  Note  that  the  color  scale  is  different 
for  each  map. 


Plate  lb.  Global  Rayleigh  wave  phase  velocity  maps  at  35,  50,  75,  and  150  s.  Note  that  the  color  scale  is 
different  for  each  map. 


Table  1.  Number  of  Observations 


Period 

N  (Love) 

N  (Rayleigh) 

35 

15,473 

28,457 

37 

15,473 

28,457 

40 

15,721 

28,663 

45 

15,780 

28,779 

50 

22,633 

37,104 

60 

23,193 

37,734 

75 

23,228 

37,739 

100 

22,498 

37,374 

150 

16,798 

33,475 
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Table  2.  Empirically  Determined  Observational  Uncertainties 


Love  Waves  Rayleigh  Waves 


Period 

CTB 

(TO 

35 

0.812 

1.164 

2.490 

1.300 

2.056 

3.519 

37 

0.727 

0.943 

1.805 

1.169 

1.753 

2.653 

40 

0.644 

0.797 

1.295 

1.045 

1.553 

2.057 

45 

0.550 

0.644 

0.922 

0.821 

1.290 

1.599 

50 

0.507 

0.709 

1.103 

0.822 

1.391 

2.016 

60 

0.388 

0.511 

0.682 

0.645 

1.042 

1.390 

75 

0.308 

0.399 

0.494 

0.444 

0.651 

0.832 

100 

0.243 

0.330 

0.410 

0.322 

0.449 

0.547 

150 

0.194 

0.253 

0.290 

0.220 

0.294 

0.310 
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Table  3.  Results  and  Comparisons  With  Other  Models 


Period 

L  = 

16^ 

L  = 

40*’ 

TW95' 

TW96'* 

M  et  al.® 

LM^ 

ZL® 

<^Crms 

Co 

7^ 

Co 

n 

^Crms 

CO 

n 

<^Crmc 

Co 

^Crms 

Co 

n 

^Crms 

Co 

n 

^Crms 

Co 

n 

L  35 

4.00 

46.2 

3.91 

78.5 

4.48 

91.0 

L  37 

3.69 

46.7 

3.62 

82.5 

4.16 

89.9 

L  40 

3.29 

47.5 

3.25 

86.9 

2.79 

105.0 

2.94 

154.2 

3.75 

88.4 

L  45 

2.83 

52.9 

2.77 

91.3 

3.21 

86.6 

L  50 

2.42 

54.3 

2.42 

101.3 

2.82 

85.3 

L  60 

2.02 

61.7 

2.02 

109.1 

2.06 

117.8 

2.05 

188.8 

2.29 

83.7 

L  75 

1.80 

70.4 

1.77 

113.0 

1.84 

82.2 

1.57 

60.1 

L  100 

1.63 

74.9 

1.56 

102.9 

1.69 

85.0 

1.72 

251.6 

1.45 

80.6 

1.41 

54.1 

1.60 

65.( 

L  150 

1.40 

72.4 

1.27 

71.7 

1.71 

73.9 

1.44 

163.3 

1.08 

78.7 

1.25 

52.1 

1.47 

55.' 

R  35 

2.49 

54.4 

2.46 

90.0 

2.00 

114.8 

R  37 

2.32 

56.5 

2.30 

92.2 

1.84 

116.8 

R  40 

2.14 

59.4 

2.13 

93.4 

2.12 

127.7 

2.08 

179.6 

1.66 

118.7 

R  45 

1.99 

62.4 

1.99 

100.2 

1.44 

119.5 

R  50 

1.90 

62.1 

1.91 

101.9 

1.30 

118.3 

R  60 

1.82 

61.7 

1.82 

101.3 

1.98 

136.9 

2.04 

200.4 

1.13 

113.2 

R  75 

1.66 

61.9 

1.69 

107.4 

1.00 

104.8 

1.43 

54.3 

R  100 

1.38 

64.1 

1.38 

103.4 

1.33 

98.1 

1.68 

247.3 

0.88 

94.3 

1.18 

48.4 

1.22 

77.; 

R  150 

0.93 

66.6 

0.87 

84.3 

0.99 

115.8 

1.10 

164.0 

0.72 

83.5 

0.92 

56.4 

0.90 

72.; 

^Undamped  inversion  with  L  =  16  (this  study). 

^Damped  inversion  with  L  =  40  (this  study). 

^Trampert  and  Woodhouse  [1995]. 

^Trampert  and  Woodhouse  [1996]. 

^Phase  velocity  maps  calculated  from  Mooney  et  al,  [1995]. 
^ Laske  and  Masters  [1996]. 

^  Zhang  and  Lay  [1996]. 
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Table  4.  Correlations  of  L  =  40  Results  With  Other  Models 


Period 

1  =  16® 
168 

TW95‘’ 

408 

TW96‘^ 

408 

M  et  al.** 
408 

LM* 

248 

IV 

308 

L  35 

0.97 

0.83 

L  40 

0.97 

0.82 

0.80 

0.79 

L  60 

0.97 

0.82 

0.78 

0.49 

I  75 

0.97 

0.28 

0.83 

L  100 

0.97 

0.83 

0.65 

0.09 

0.81 

0.82 

L  150 

0.95 

0.81 

0.67 

-0.09 

0.85 

0.82 

R35 

0.98 

0.54 

R  40 

0.98 

0.84 

0.81 

0.36 

R  60 

0.98 

0.84 

0.79 

-0.10 

R  75 

0.98 

-0.22 

0.81 

R  100 

0.99 

0.80 

0.68 

-0.28 

0.81 

O.OO'’ 

R  150 

0.98 

0.64 

0.62 

-0.23 

0.77 

0.62‘ 

*  Undamped  inversion  with  L  =  16  (this  study). 

^Ttamperi  and  Woodhouse  [1995]. 

^Trampert  and  Woodhouse  [1996]. 

^  Phase  velocity  maps  calculated  from  Mooney  et  al.  [1995]. 
^Laske  and  Masters  [1996]. 

^  Zhang  and  Lay  [1996]. 

^Highest  degree  (Lm&x)  in  model. 

=  28. 

■Ln,ax  =  24. 
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